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ABSTRACT 


Until  recently,  computer  simulations  of  helicopter  rotor  dynamics  have  employed 
equations  of  motion  that  have  been  linearized  or  simplified.  These  modified  equations  of 
motion  did  not  allow  for  the  evaluation  of  nonlinear  material  properties  in  the  rotor  since 
higher  order  terms  in  the  dynamics  had  been  modified  in  the  simplification  process.  With 
recent  advances  in  both  computer  simulation  hardware  and  symbolic  mathematic 
manipulation  software,  the  full  nonlinear  equations  of  motion  may  be  utilized  in  helicopter 
rotor  simulations.  This  dissertation  reports  on  the  use  of  the  full  nonlinear  equations  of 
motion  in  the  analysis  of  rotor  blade  lead/lag  motion  and  its  effect  on  rotor  hub  and  rigid 
body  fuselage  motion.  Nonlinear  modeling  methods  are  implemented  using  Maple 
symbolic  mathematic  manipulation  software  and  Matlab  and  Simulink  computer 
simulation  environments.  Results  are  compared  to  the  RAH-66  Comanche  Froude  scale 
wind  tunnel  article  and  new  methodologies  evaluated  in  the  search  for  a  damperless  rotor 
system  that  is  free  of  ground  and  air  resonance  mechanical  instabilities. 
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I.  INTRODUCTION 

One  of  the  goals  of  this  research  was  to  develop  a 
flexible  computational  tool  to  analyze  the  dynamic  and 
aeromechanical  behavior  of  advanced  technology  coupled 
rotor/fuselage  systems.  Initially,  a  series  of  programs  were 
developed  utilizing  the  symbolic  processing  software. 
Maple®,  the  computational  software,  Matlab®,  and  the 

simulation  software,  Simulink®,  by  LT  Christopher  Robinson 
[Ref.  1].  It  was  desired  that  the  computational  tool  be 
simple  to  understand  and  lend  itself  to  easy  reprogramming 
by  any  user  knowledgeable  in  the  field  of  dynamics  and  in 
the  use  of  the  software  programs  mentioned.  It  was  also 
desired  that  the  developed  programs  allow  for  a  sufficient 
capability  so  that  the  effects  of  introducing  advanced 
technologies  into  rotor  system  designs  could  be  accurately 
modeled.  This  dissertation  reports  on  further  advances  in 
this  computational  tool  development,  and  analysis  based  on 
these  new  tools.  The  nonlinear  rotor  simulation  that  was 
developed  as  part  of  this  research  is  a  very  powerful  tool 
for  analyzing  rotors  that  do  not  have  multiple  load  paths  to 
the  hub. 

Historically,  rotor  analysis  was  performed  using 
approximate  equations  of  motion.  In  order  to  reduce  the 
complexity  of  the  computer  code  involved  and  to  reduce 
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computation  times,  these  equations  were  simplified  by 
eliminating  higher  order  terms  [Refs.  2,3].  Subsequently, 
either  linearized  equations  were  used  or  select  nonlinear 
terms  were  retained  using  a  ranking  system  or  ordering 
scheme,  such  as  that  used  by  Friedmann,  [Ref.  4]. 

As  hingeless  helicopter  main  rotors  became  more 
commonplace,  a  need  arose  for  a  rotor  simulation  tool  that 
would  accurately  model  nonlinear  mechanical  properties  so 
that  these  nonlinearities  could  be  exploited.  It  was  desired 
to  develop  an  analysis  tool  that  would  be  able  to  reliably 
model  the  effects  of  nonlinearities  in  the  rotor  and  hub 
mechanical  and  geometric  parameters.  For  this  reason, 
Robinson  and  Wood  developed  a  utility  for  creating  the  full 
nonlinear,  coupled  equations  of  motion.  Robinson  utilized 
the  symbolic  manipulation  software  Maple®  for  this  purpose. 

Current  trends  in  helicopter  technology  and 
manufacturing  have  favored  the  use  of  bearingless  rotor 
designs  that  make  use  of  advanced  composite  materials.  These 
designs  offer  many  advantages  over  more  conventional 
articulated  rotors  in  reliability  and  maintainability.  A 
potential  payoff  from  the  successful  use  of  the  technologies 
mentioned  above  is  the  damperless  rotor;  a  design  that 
offers  major  returns  in  the  form  of  decreased  rotor  system 
weight,  reduced  parts  count,  and  reduced  maintenance 
requirements . 
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Composites  and  other  advanced  materials  that  can  be 
applied  for  a  damperless  rotor  make  modeling  rotor  behavior 
more  difficult,  however.  These  materials  have  the  potential 
for  exhibiting  nonlinear  behavior  that  cannot  be  accurately 
modeled  without  utilizing  the  full  nonlinear  equations  of 
motion  for  the  system. 

The  goal  of  this  dissertation  is  to  examine  alternative 
designs  for  a  damperless  helicopter  rotor.  This  effort  met 
with  success  by  taking  advantage  of  nonlinear  stiffness  in 
the  blade  root  end  to  avoid  divergent  motion  in  the  rotor 
blade  lead/ lag  degrees  of  freedom.  This  oscillatory 
instability  is  better  known  as  ground  and  air  resonance 
instability. 
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II.  REVIEW  OF  WORK  PERFORMED  BY  ROBINSON 

For  background  on  progress  reported  in  this 
dissertation,  it  is  important  to  first  review  LT  Robinson's 
thesis  work.  This  work  preceded  that  of  the  author.  Tasks 
performed  by  Robinson  included  the  formulation  of  a  Maple® 

based  symbolic  processing  worksheet  that  formulated 
nonlinear  equations  of  motion  given  energy  expressions  for 
helicopter  rotor  model  degrees  of  freedom.  Simulink®  based 
computer  simulations  were  developed  from  the  equations  of 
motion  derived  by  the  symbolic  processor  for  a  simple  three 
bladed  rotor  based  on  that  used  by  Coleman  [Refs.  5,6] . 

Robinson,  Wood,  and  King  reported  on  a  new  method  for 
formulating  the  full  non-linear  equations  of  motion  for 
ground/air  resonance  stability  analysis  of  helicopter  rotor 
systems  [Ref.  7],  The  full  set  of  non-linear  equations  was 
developed  by  Lagrangian  approach  using  symbolic  processing 
software  for  expanding  the  equations.  The  symbolic  software 
was  further  utilized  to  automatically  convert  the  equations 
of  motion  into  C  or  Fortran  source  code  formatted  for 
numerical  integration.  Simulink®  then  applied  a  Runga-Kutta 

integration  scheme  to  generate  time  history  plots  of  blade 
and  fuselage  motion.  Damping  levels  were  determined  from  the 
time  history  simulations  by  a  Matlab®  program,  which  used 
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the  Moving  Block  Technique  for  establishing  damping  levels 
from  the  coupled  rotor-body  response. 

Robinson's  model,  based  on  that  of  Coleman,  is  shown  in 
Figure  1,  [Ref.  1] . 


Figure  1  -  Simplified  Rotor  Model  [After  Ref.  1] 

Robinson  verified  his  symbolic  Lagrangian  derivation  by 
comparing  his  results  to  Coleman's.  Figure  2  shows  results 
of  a  comparison  of  the  two  methods  for  an  unstable  case 
where  a  moderate  amount  of  damping  has  been  added  to  both 
rotor  blades  and  fuselage. 
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time  (sec) 


Figure  2  -  Validation  of  simple  model  with  solution  of  Coleman’s  linearized  model 


In  Figure  2,  Robinson  shows  excellent  agreement  between 
the  two  solutions  with  departure  occurring  only  when 
displacements  are  very  large.  This  is  to  be  expected  since 
the  simulation  model  of  this  paper  does  not  assume  small 
angle  theory  whereas  the  Coleman-Feingold-Bramwell  model 
does . 

For  a  more  quantitative  proof  of  the  validity  of  the 
new  method,  it  was  tested  against  Deutsch's  Criteria  [Ref. 
8].  Based  on  Coleman's  analysis,  Deutsch  showed  that  the 
product  of  the  hub  (landing  gear)  damping  and  the  blade 
lead/ lag  damping  must  be  greater  than  a  prescribed  value  to 
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collapse  the  band  of  unstable  blade  frequencies :  A,p  Xj,  >  A,3 
(p-1) .  The  following  Figure  shows  results  of  successive 
sweeping  across  the  unstable  band  using  the  new  nonlinear 
analysis.  The  respective  cases  represent  those  where 
Deutsch's  prescribed  damping  is  applied  (D=1.0),  and  other 
cases  where  the  damping  is  greater  than  that  required  by 
Deutsch  (D=l .  1 )  or  less  (D=0.8,  0.9).  For  each  case 
considered,  the  center  of  instability  is  seen  to  be  the 
point  of  minimum  value  on  the  given  line.  With  the  new 
nonlinear  analysis,  Deutsch's  criteria  applied  at  the  center 
of  instability  provides  nearly  a  neutrally  stable  result 
with  a  near  zero  damping  ratio.  Deutsch's  criteria  is 
conservative  in  that  the  critical  point  for  his  criteria,  or 
bucket  of  each  curve,  still  shows  positive  damping  when 
analyzed  with  the  NPS  simulation.  Three  years  after 
Robinson,  Rafanello  verified  the  conservative  nature  of 
Deutsch's  criteria  by  matching  simulation  results  to  Wood's 
HSS-2  ground  resonance  analysis,  [Refs.  9,  10] . 
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Figure  3  -  Quantitative  validation  of  simple  model  against  Deutsch  Criteria 
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III.  DAMPERLESS  ROTOR  DEVELOPMENT 


A.  BACKGROUND 

Ground  and  air  resonance  are  particularly  destructive 
mechanical  instabilities  that  can  occur  in  helicopters  with 
fully  articulated,  bearingless,  or  hingeless  main  rotor 
designs.  The  phenomenon  of  ground/air  resonance  is  the 
result  of  coupling  between  rigid  body  motions  of  the 
fuselage  on  its  landing  gear,  or  in  the  air,  and  lead-lag 
oscillations  of  the  rotor  blades  in  their  plane  of  rotation. 
When  it  occurs,  the  instability  can  build  to  destructive 
proportions  in  a  matter  of  seconds . 

The  equations  of  motion  describing  the  dynamics  of  the 
coupled  rotor-fuselage  system  are  nonlinear  and  generally 
quite  complex  even  for  simplified  models.  The  fundamental 
theory  of  air  and  ground  resonance  is  credited  to  the 
classic  work  of  an  NACA  engineer,  Robert  Coleman,  who  first 
published  his  theory  in  the  early  1940's  [Ref.  5] .  As 
computational  power  improved  with  the  evolution  of  digital 
computers,  more  general  techniques  for  analyzing  rotor 
system  stability  were  developed. 

As  stated  earlier,  due  to  the  complexity  of  the 
equations  of  motion,  ground  and  air  resonance  analysis  was 
historically  restricted  to  linear  analysis  until  very 
recently.  This  frequently  restricted  the  scope  of  many 
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studies  and  potential  solutions.  With  a  few  noteworthy 
exceptions,  nonlinearities  in  ground  and  air  resonance 
research  programs  have  not  typically  investigated  beyond 
traditional  damping  methods,  such  as  hydraulic  damping. 

This  chapter  entitled  "Damper less  Rotor  Development" 
investigates  various  options  for  damperless  rotors  that 
maintain  some  margin  of  stability,  or  avoidance  of  ground 
and  air  resonance  without  the  employment  of  auxiliary 
lead/ lag  dampers. 

In  general,  rotorcraft  lag  dampers  are  heavy,  expensive 
items  that  can  require  extensive  maintenance.  By  developing 
rotors  that  do  not  require  these  dampers,  payoffs  in  1) 
reduced  maintenance  time  and  cost,  2)  reduced  hub  complexity 
and  parts  count,  3)  reduced  rotor  drag,  and  4)  reduced  rotor 
signature  could  potentially  be  reaped. 

An  effort  in  the  elimination  of  these  dampers  applies 
directly  to  a  famous  quote  from  the  father  of  the  business 
jet.  Bill  Lear,  "Strive  for  design  simplicity:  you  never 
have  to  fix  anything  you  leave  out." 


B.  GROUND  RESONANCE  FUNDAMENTALS 

Ground  and  air  resonance  instability  is  centered  around 
blade/hub  interactions  at  one  of  the  fuselage  rigid  body 
modes.  The  rotor  acts  as  a  forcing  function,  producing 
undesired  responses  or  even  resonance  in  the  fuselage 
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dynamics .  The  rotor  changes  the  frequency  of  mechanical 
vibrations  as  they  transfer  from  the  rotating  frame  down 
into  the  non-rotating  frame.  If  the  blades  have  a  primary 
mode  in  lead/lag  at  some  frequency  (0lag,  and  £2  is  the 
frequency  of  rotation  of  the  rotor  system,  shears  are 
transmitted  through  the  mast  and  into  the  non-rotating  hub 
at  frequencies  (£2±colag)  .  The  progressing  mode  (£2+(0lag)  is 
heavily  damped  and  normally  plays  no  role  in  ground 
resonance.  The  regressing  mode  (£2-CDlag)  provides  the 
potential  for  destruction,  however.  If  one  of  the  fuselage 
rigid  body  modes  is  near  in  frequency  to  the  regressing  mode 
and  there  is  insufficient  damping,  blade  lead/ lag  and 
fuselage  displacement  amplitudes  can  build  to  destructive 
proportions.  While  the  helicopter  is  still  in  contact  with 
the  ground,  the  fuselage  roll  mode  is  often  close  in 
frequency  to  the  regressing  lag  mode,  normally  in  the  2-5  Hz 
range . 

In  general,  the  spectrum  of  classical  ground  resonance 
to  air  resonance  could  be  referred  to  as  " aeromechanical 
stability."  Since  articulated  rotors  normally  do  not 
experience  air  resonance,  aerodynamics  and  flap  degrees  of 
freedom  play  little  roll  in  ground  resonance.  Therefore, 
articulated  rotors  may  be  modeled  using  Coleman's  methods. 
Hingeless  and  bearingless  rotors  experience  more  coupling  in 
the  blade  degrees  of  freedom.  In  addition,  body  pitch  and 
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roll  coupling  through  flap  moments  of  cantilever  blades  can 
result  in  air  resonance.  Bearingless  and  hingeless  rotors, 
where  Coleman's  equations  are  not  sufficient  to  model  the 
problem,  require  more  sophisticated  models.  Flap  degrees  of 
freedom  and  aerodynamics  were  not  used  to  model  the  Coleman¬ 
like  rotor  models  in  this  report. 

In  the  typical  rotor,  there  are,  in  fact,  significant 
nonlinearities.  Especially  noteworthy  are  friction  in 
dampers  and  oleos  that  sometimes  result  in  the  appearance  of 
instability  following  a  moderately  large  external  force 
excitation.  Thus,  the  typical  scenario  for  ground  resonance 
is  initiated  by  the  fuselage  being  perturbed  by  some 
external  force  such  as  a  gust,  a  rolling  ship  deck,  uneven 
terrain  while  ground  taxing,  or  combat  damage.  The  fuselage 
motion,  typically  in  roll,  induces  asymmetrical  lag  motion 
in  the  rotor  blades.  This  asymmetrical  response  shifts  the 
combined  CG  of  the  rotating  components  away  from  the  center 
of  rotation,  resulting  in  an  unbalance  in  the  centrifugal 
force  seen  by  the  rotor  mast.  The  unbalanced  centrifugal 
force  results  in  increased  fuselage  motion,  compounding  the 
problem.  This  cycle  is  depicted  in  flowchart  form  in  Figure 
4. 
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Figure  4  -  Ground  resonance  flowchart 


A  typical  fuselage  on  its  landing  gear  has  a  roll  mode 
frequency  in  the  proximity  of  the  regressing  mode  frequency 
of  the  rotor  system.  If  these  modes  are  sufficiently  near  to 
each  other  and  there  is  insufficient  system  damping,  ground 
resonance  can  occur. 
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A  plot  of  hub  motion  for  a  typical  case  of  ground 
resonance  can  be  found  in  Figure  5.  Note  the  spiraling 
divergence  of  the  hub  center  of  mass  in  the  horizontal 
plane . 


-1  -0.8  -0.6  -0.4  -0.2  0  0.2  0.4  0.6  0.8 

Hub  X  displ  (ft) 

Figure  5  -  Ground  resonance  hub  motion 

Following  the  work  of  Robert  Coleman,  a  more  complete 
work,  based  on  his  analysis,  was  collected  and  published  by 
Coleman  and  Feingold  [Ref.  6].  Deutsch  supplemented 
Coleman's  work  by  establishing  a  criterion  for  system 
damping  necessary  to  eliminate  ground  resonance  [Ref.  8] . 
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Numerous  other  researchers  have  contributed  to 
understanding  the  ground  resonance  problem  since  Coleman, 
Deutsch,  and  Feingold.  The  papers  following  Coleman 
essentially  completed  the  investigation  of  the  linearized 
system  for  ground  resonance.  Hammond  applied  Floquet  theory, 
but  to  consider  the  case  of  unequal  lead/ lag  damping  among 
the  blades  [Ref.  11]  .  (The  NPS  rotor  simulation  can  also 
apply  unequal  lead/lag  damping  values  at  each  blade.) 
Tongue,  Flowers,  Jankowski,  Tang,  Dowell,  and  most  recently, 
Gandhi  and  Chopra  [Refs.  12-14]  investigated  nonlinear 
effects  in  blade  and  fuselage  damping.  Ormiston  explored 
linear  aeromechanical  stability  of  rigid  blades  with  spring 
restrained  hinges,  and  later,  investigated  rotor  modeling 
with  a  sophisticated  finite  element  analysis  code  and 
various  nonlinear  damping  models  [Refs.  15,  16]. 


C.  OBJECTIVES 

The  objective  of  the  present  work  is  to  explore  the 
potential  of  eliminating  the  snubber-damper  or  damper  on 
hingeless  rotor  designs  and  replace  it  with  a  f lexbeam  that 
has  been  modified  to  possess  nonlinear  properties.  Thus,  the 
purpose  is  to  employ  the  nonlinear  properties  of  the 
flexbeam  to  introduce  nonlinear  dynamic  characteristics  that 
will  replace  unbounded  instability  with  small  amplitude 
limit  cycle  oscillations,  and  prevent  catastrophic 
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destruction  of  the  helicopter.  This  work  is  well  suited  to 
the  new  NPS  rotor  simulation  analysis  that  can  accurately 
model  nonlinear  mechanical  properties. 

The  significance  of  this  new  research  is  that  allowance 
of  nonlinearities  at  the  blade  root  may  result  in  an 
acceptable  bounded  response  in  the  parameter  region  where 
linear  theory  would  predict  instability,  or  may  introduce 
new  regions  of  unacceptable  response  in  the  parameter  region 
where  linear  theory  would  predict  stability.  Evidence  of  the 
latter  may  be  found  in  the  numerous  aircraft  lost  in  the 
documented  cases  of  ground  resonance.  Modern  soft-inplane 
rotors,  such  as  that  first  introduced  on  the  BO-105,  have 
the  additional  possibility  of  encountering  this  lead  lag 
instability  in  flight  [Ref.  17] . 

The  work  reported  here  attempts  to  abandon  the  use  of 
auxiliary  damping  methods  in  the  rotating  system  as  a  means 
of  eliminating  ground  and  air  resonance.  Only  the  damping 
inherent  to  the  typical  landing  gear  system  and  the  built-in 
structural  damping  of  the  flexbeam  in  bending  are  utilized 
in  the  cases  shown.  The  ground/air  resonance  cycle  is 
detuned  by  shifting  the  blade  natural  frequency  in  lead/lag 
to  avoid  the  coalescence  of  the  regressing  mode  with  the 
fuselage  rigid  body  modes .  A  nonlinear  flexbeam  is 
characterized  by  change  in  stiffness  properties  with  change 
in  amplitude.  Therefore,  as  blade  lead/lag  amplitude 
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changes,  there  is  a  corresponding  change  in  its  lead/lag 
frequency.  This  method  differs  from  some  previous  methods  of 
preventing  ground  resonance  in  that  coalescence  is  avoided 
through  natural  frequency  shift  instead  of  dominating  the 
dynamics  of  the  system  through  damping.  In  the  case  of  this 
thesis,  this  is  accomplished  by  introducing  a  classic 
nonlinear  stiffness  known  as  Duffing-type  nonlinear 
stiffness  at  the  blade  root. 

D.  INTRODUCTION  TO  DUFFING 

1.  Single  Degree  of  Freedom  Duffing  Systems 

For  the  first  part  of  the  investigation,  an  exploration 
of  the  one  degree  of  freedom  Duffing  system  is  considered  to 
understand  the  effects  of  cubic  stiffness  for  a  simple 
harmonic  oscillator.  With  damping  excluded.  Duffing' s 
equation  can  be  expressed  as : 

x  +  co^x  +  hx^  =Gcos(fi)jO 

where  the  driving  function  is  Gcos^t)  .  The  coefficient,  G, 
is  the  ratio  between  the  amplitude  of  the  force  and  the  mass 
of  the  moving  system.  Coefficients  co02  and  h  depend  on  the 
parameters  of  the  system,  with  the  Duffing  coefficient,  h, 
being  positive  for  a  "hardening  spring"  and  negative  for  a 
"softening  spring." 
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One  of  the  interesting  phenomena  appearing  in  systems 
governed  by  Duffing' s  equation  is  that,  in  addition  to 
generating  harmonics  of  the  driving  frequency,  discontinuous 
jumps  in  amplitude  occur  as  the  frequency  of  the  exciting 
force  is  either  increased  of  decreased  at  a  constant  rate. 

A  simple  one  degree  of  freedom  Duffing  spring  mass 
system  is  shown  in  Figure  6. 


x(t) 


F  *  sin(w*t) 

Figure  6  -  Simple  Duffing  system 

A  simple  simulation  was  written  to  obtain  time 
histories  for  this  system  with  user  supplied  mass,  linear 
damping,  linear  stiffness,  forcing  amplitude,  forcing 
frequency,  and  Duffing  stiffness  terms.  The  Simulink®  5th 
order  Runge-Kutta  integration  scheme  is  used  here. 


20 


An  example  time  history  that  illustrates  the  effect  of 
exciting  frequency  (both  increase  and  decrease)  is  presented 
in  Figure  7 . 


Duffing  System  Time  Histories 
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Figure  7  -  Duffing  time  history 

The  inputs  for  this  time  history  are  contained  in  the 
following  table: 


Table  1  •  Table  caption 


Linear  stiffness  =  4000  lb/in 

Mass  =  6  slugs 

Duffing  stiffness  =  100  lb/inJ 

Damping  =  14  lb*s/in 
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Driving  Force  1/2  amplitude  =  500  lb 


The  top  subplot  in  Figure  7  is  the  time  history  of  the 
mass  displacement  with  a  stiffening  Duffing  spring  included 
in  the  dynamics .  The  second  subplot  in  this  figure  is  the 
driving  frequency  used  in  the  top  plot  as  a  function  of 
time.  The  system  goes  through  resonance  both  as  the  driving 
frequency  is  increased  and  again  as  it  decreases. 

If  one  takes  this  time  history  and  transforms  it  as  a 
function  of  driving  frequency,  instead  of  time, 
dissimilarities  in  the  two  resonances  can  easily  be 
discerned,  as  shown  in  Figure  8. 
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Mass  Deflection  Time  Histories 


Frequency  (Hz) 


Figure  8  -  Duffing  system  resonance 

The  system  reaches  resonance  as  the  frequency  is 
increased,  but  unlike  a  linear  system  that  gradually  reduces 
the  displacement  amplitude,  in  this  case,  the  Duffing  system 
experiences  a  sudden  drop-off  in  displacement  amplitude. 
Then,  as  the  driving  frequency  is  decreased  from  a  value 
above  resonance,  the  system  does  not  follow  the  same  path  in 
resonance  amplitude,  but  reaches  a  relative  maximum  in 
displacement  at  a  frequency  less  than  the  case  for 
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increasing  driving  frequency.  These  sudden  changes  in 
displacement  amplitude  are  referred  to  as  "Duffing  jumps." 

Shown  in  Figure  9  is  a  nonlinear  resonance  curve  for  a 
hardening  Duffing  spring  system,  reproduced  from  Cunningham 
[Ref.  18] .  A  typical  linear  resonance  curve  has  been 
superimposed  on  Cunningham's  original  figure.  The  resonant 
frequency  for  the  linear  system  would  not  change,  but  remain 

fixed  at  cd0.  The  physical  reason  for  these  jumps  can  be 

visualized  in  this  plot  of  the  resonance  of  the  nonlinear 
Duffing  system. 

As  the  driving  frequency  is  increased  from  the  left  on 
the  plot,  the  system  displacement  amplitude  grows  until  the 
system  reaches  point  A.  Further  increase  in  the  driving 
frequency  drops  the  displacement  amplitude  down  to  B,  a 
Duffing  jump. 

As  the  driving  frequency  is  decreased  from  the  right, 
the  system  drives  right  through  point  B  to  point  C,  where 
the  solution  becomes  unstable.  As  the  driving  frequency  is 
decreased  lower  than  that  found  at  point  C,  the  only  stable 
displacement  amplitude  becomes  point  D,  another  Duffing 
jump. 
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Figure  9  -  Duffing  resonance  frequencies  [After  Ref.  18] 


Unlike  the  linear  system,  the  Duffing  system  increases 
the  natural  frequency  as  a  function  of  displacement 
amplitude.  For  larger  displacements,  the  resonance  frequency 
increases  due  to  this  stiffening.  The  relation  for  this 
frequency  shift  is  derived  by  Cunningham  as: 


It  should  be  noted  that  if  the  sign  of  h  is  allowed  to 
reverse,  h<l,  (softening  spring),  the  same  kind  of  phenomena 
exists.  The  primary  difference  is  that  now  the  nonlinear 


resonance  curve  is  skewed  towards  the  lower,  rather  than 
higher  frequencies .  Since  ground  and  air  resonance  are 
determined  by  coupling  of  the  rotor  blades  with  lead/ lag 
motion  in  the  plane  of  rotation,  it  was  hoped  to  take 
advantage  of  these  sudden  decreases  in  displacement 
amplitude.  This  decrease  in  blade  lead/lag  motion  would 
subsequently  reduce  the  severity  of  the  offset  centrifugal 
force  on  the  hub.  The  employment  of  Duffing  in  the  rotor 
lead/lag  dynamics  did,  in  fact,  reduce  the  susceptability  of 
the  rotor  to  ground  and  air  resonance,  but  not  due  to  the 
Duffing  jumps  as  hypothesized.  The  shift  in  natural 
frequency  detuned  the  ground  resonance  sequence  of  events, 
instead. 


2.  Multiple  Degree  of  Freedom  Duffing  Systems 

Since  the  rotor  system  to  be  simulated  was  a  three 
bladed  design,  investigations  in  a  three  degree  of  freedom 
Duffing  system  were  also  performed.  A  depiction  of  this 
system  may  be  found  in  Figure  10 .  The  possibility  existed  in 
this  multiple  DOF  nonlinear  system  that  reductions  in 
resonant  displacement  amplitude  could  be  achieved  through 
coupling  between  the  three  degrees  of  freedom  in  the 
helicopter  rotating  frame. 
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Figure  10  -  Three  DOF  Duffing  system 

As  shown,  the  Duffing  spring/mass  system  described 
earlier  was  expanded  into  a  three  mass  system.  Each  mass  was 
linked  by  adjacent  linear  springs,  linear  dampers,  and 
Duffing  springs.  An  example  time  history  of  this  three  DOF 
system  follows.  The  input  parameters  are  the  same  as  the  one 
DOF  case  shown  previously.  The  three  resonant  peaks  are 
clearly  visible. 
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Figure  11  -  Three  DOF  Duffing  system  time  histories 


In  general,  no  substantial  changes  in  the  dynamics  in 
the  neighborhood  of  the  first  resonant  frequency  of  this 
system  were  noted  over  the  one  DOF  system.  The  second  and 
third  modes  were  most  effected  by  variation  in  the  Duffing 
terms  of  the  system.  By  increasing  the  Duffing  terms,  the 
higher  order  modes  tended  to  flatten  out,  showing  resonant 
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behavior  over  larger  ranges  of  the  driving  frequency.  Three 
plots  showing  this  trend  are  shown  in  Figure  12,  first  the 
linear  system  with  no  Duffing  stiffness,  with  an 
intermediate  Duffing  coefficient,  then  with  a  large  Duffing 
coefficient. 


Increasing  Dulling  Stiffness  in  3  DOF  Duffing  System 


Time  (sec) 


0  2  4  6  8  10  12  14  16  18 

Freq  (Hz) 

Figure  12  -  Effect  of  increasing  Duffing  stiffness  in  three  DOF  system 
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Since  the  ground  and  air  resonance  problem  is  more 
closely  associated  with  the  lower  resonant  frequencies  of 
the  rotor,  no  advantage  in  the  ground/air  resonance  problem 
was  gleaned  from  these  higher  order  Duffing  system  analyses. 
The  shift  in  resonance  frequencies,  however,  still  showed 
promise  in  rotor  applications,  as  discussed  in  the  following 
section. 


E.  DUFFING  ROTOR  SYSTEMS 

As  discussed  previously,  ground  and  air  resonance  are 
effectively  driven  by  two  phenomena:  1)  the  proximity  of 
rotor  and  hub  vibratory  modes,  and  2)  the  ability  of  system 
damping  to  maintain  reasonable  DOF  amplitudes  when 
approaching  resonance.  Damping,  the  traditional  means  of 
solving  the  ground  resonance  problem,  can  dominate  the 
system  dynamics  through  the  velocity  terms.  Except  for  high 
damping  ratios,  this  approach  does  little  to  shift  system 
modal  frequencies ,  however . 

The  method  of  passive  rotor  stability  presented  here 
utilizes  a  nonlinear  Duffing  type  spring  to  address  the 
coalescence  part  of  the  problem,  while  attempting  to  settle 
for  whatever  damping  is  inherent  to  the  system.  This 
approach  eliminates  the  need  for  auxiliary  damping  devices 
for  the  blade  lag  degrees  of  freedom. 
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Duffing  stiffness  is  added  into  the  Coleman  model  as 
depicted  in  Figure  13 .  Both  linear  and  nonlinear  rotary 
springs  plus  structural  damping  are  incorporated  across  the 
blade  lag  hinge.  Stiffness  and  damping  values  can  be  defined 
separately  for  each  blade  to  allow  for  asymmetric  blade 
property  investigations . 

The  addition  of  a  hardening  Duffing  spring  increases 
the  lead/ lag  natural  frequency  of  the  blade  as  a  function  of 
blade  lead  angle,  detuning  the  resonance.  As  blade 
amplitudes  increase  in  a  typical  ground  resonance  sequence 
of  events,  the  regressing  mode  is  driven  to  a  lower 
frequency  as  the  blade  natural  frequency  in  lag  is  increased 
due  to  stiffening.  This  effectively  breaks  the  sequence  of 
events  previously  depicted  in  flowchart  form.  This  shift  in 
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the  rotor's  regressing  mode  puts  the  blade  lag  motion  in  a 
limit  cycle  instead  of  letting  displacements  increase 
without  bound.  For  the  case  of  soft  in-plane  rotors  that 
drive  through  the  ground  resonance  region  of  instability 
upon  rotor  engagement,  upon  leaving  this  region,  the 
dynamics  return  to  what  one  would  expect  from  a  soft  in¬ 
plane  rotor  at  rotor  speeds  above  the  region  of  instability. 
Outside  the  region  of  instability,  the  linear  terms  once 
again  dominate  the  dynamics,  with  little  or  no  ill  effects 
of  the  Duffing  stiffness  seen  in  the  dynamics . 

1.  Results  of  the  Baseline  Duffing  Rotor 

The  NPS  nonlinear  rotor  simulation  was  utilized  in  the 
following  investigation.  To  provide  the  most  robust 
application  of  Duffing  springs  to  the  ground  resonance 
problem,  simulations  are  performed  at  the  center  of  the 
region  of  instability  unless  otherwise  stated.  This  baseline 
case  is  of  a  rotor  with  the  following  values: 

Table  2  -  Table  caption 


Q  =  25.0  rad/s 
(239  RPM) 

Mhub  =  180.0  si 
(5800  lbs) 

M.i.ae  =  6.0  si 
(193  lbs) 

e  =  0.5  ft 

Rbiadecc  =  10-5  ft 

(0lag/  Q  =  0.6 

Structural  Damping  in  Lag  =  2.0% 

Viscous  Damping  in  Hub  =2.0% 
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Time  histories  of  blade  and  hub  displacements  for  this 
rotor  with  no  Duffing  stiffness  employed  are  illustrated  in 
Figure  14  and  Figure  15 .  The  first  plot  shows  the  blade  lag 
angles  clearly  diverging  from  their  initial  zero 
displacement  angles . 


012345678 


Time  (sec) 

Figure  14  -  Baseline  case,  blade  lead  angles  without  Duffing  stiffness 

The  X  direction  of  the  isotropic  hub  motion  is  plotted 
in  Figure  15.  The  hub  has  an  initial  lateral  displacement  of 
0.2  feet  to  initiate  the  ground  resonance  behavior  of  the 
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system.  Like  the  blade  plot  above,  the  hub  motion  clearly  is 
subject  to  an  oscillatory  divergence. 


2.  Baseline  Case  with  Duffing  Stiffness  Added 

The  next  set  of  time  histories  gives  results  for  the 
same  rotor,  but  with  Duffing  spring  terms  added.  The  ratio 
of  the  Duffing  stiffness  coefficient  to  the  linear  stiffness 
coefficient  is  20.0  for  Figure  16  and  Figure  17. 
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Figure  16  -  Baseline  case,  blade  lead  angles  with  Duffing  stiffness 

The  nonlinear  stiffness  clearly  keeps  the  previously 
unstable  blade  lag  response  in  check.  Since  the  rotor  for 
this  case  is  linearly  unstable,  a  limit  cycle  of 
approximately  8  degrees  in  the  blade  lag  motion  is  produced. 

A  plot  of  the  hub  motion  follows.  Limit  cycle  motion, 
in  general,  is  not  considered  desirable  because  it  results 
in  additional  structural  vibration  in  the  non-rotating 
aircraft  components.  The  amplitude  of  the  lateral  motion 
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depicted  in  this  plot  is  sufficiently  large  that  concerns 
over  hub  vibration  levels  were  raised. 


Time  (sec) 


Figure  17  -  Baseline  case,  hub  lateral  displacement  with  Duffing  stiffness 

The  lateral  vibrations  experienced  by  the  hub  are 
plotted  in  Figure  18.  The  derivative  of  the  hub  velocity  was 
taken  numerically  to  obtain  the  hub  lateral  vibration  data 
plotted.  Since  numerical  derivation  of  a  sampled  signal 
tends  to  increase  the  effect  of  noise  in  the  data,  this  plot 
is  provided  as  a  means  of  discerning  the  overall  vibration 
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levels,  not  as  an  explicit  means  of  calculating  the 
vibration  time  history.  The  final  vibratory  levels  are  on 
the  order  of  one-half  the  acceleration  of  gravity. 
Vibrations  of  this  magnitude  might  be  larger  than  desired, 
even  for  transient  ground  operations. 


Time  (sec) 

Figure  18  -  Hub  lateral  vibration,  baseline  case  with  Duffing  stiffness 
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3.  Baseline  Case  with  Increased  Rotor  Speed  -  Above 
the  Region  of  Instability 

The  case  that  produced  the  unacceptable  vibration 
levels  presented  above  is  one  that  was  originally  unstable 
and  represents  a  worst-case  operating  regime  for  the  rotor: 
the  center  of  the  region  of  instability.  Even  if  the 
nonlinear  stiffness  was  included  in  the  rotor  design,  it  is 
doubtful  a  production  rotor  would  be  designed  to  operate  at 
its  center  of  instability. 

Most  soft  in-plane  rotors  are  designed  to  operate  with 
rotor  speeds  slightly  above  the  region  of  instability.  The 
rotor  drives  through  the  unstable  region  on  engagement  and 
disengagement.  It  is  also  possible  to  encounter  these  low  Q 

regimes  during  extreme  maneuver,  maintenance  checks  or  as 
the  result  of  a  system  failure. 

For  the  baseline  rotor  presented  here,  the  blade  center 
of  mass  is  at  a  radial  location  of  10.5  feet;  it  roughly 
represents  a  20-ft  radius  rotor.  For  this  size  blade,  25 
rad/sec  (239  RPM)  is  an  unusually  low  operating  speed.  If  £2 
is  increased  to  approximate  a  650  ft/sec  tip  speed  in  a 
hover,  Q  would  increase  to  32.5  rad/sec  (310  RPM). 

To  examine  the  effects  of  the  inclusion  of  nonlinear 
stiffness  at  the  higher  rotor  speed,  plots  of  the  blade  and 
hub  motion  are  provided  in  Figure  19  and  Figure  20  for 
comparison  with  the  linear  rotor  shown  previously. 
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Figure  19  -  Blade  lead  angles,  increased  rotor  speed  with  Duffing  stiffness 

Since  the  rotor  is  now  being  operated  with  a  rotor 
speed  slightly  above  the  region  of  instability,  the  hub  and 
blade  displacements  are  less  for  the  same  initial 
conditions.  The  lag  limit  cycle  amplitude  is  reduced  from  8 
degrees  all  the  way  to  0.3’ degrees. 
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Base  Case  w/  Duffing  &  Incr  Omega 


Time  (sec) 


Figure  20  -  Hub  lateral  displacement,  increased  rotor  speed  with  Duffing  stiffness 

The  hub  lateral  displacement  saw  a  substantial 
decrease,  as  well.  The  hub  limit  cycle  amplitude  for  the 
Duffing  rotor  drops  80%  from  the  value  it  had  at  the  lower 
rotor  speed. 

For  comparison,  the  hub  lateral  vibrations  for  the 
Duffing  rotor  at  both  rotor  speeds  are  plotted  in  Figure  21. 
The  vibration  amplitude  drops  almost  an  order  of  magnitude 
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Hub  Vibrations  (g’s) 


to  less  than  0.05  g's,  a  level  that  is  certainly  in  the 
window  of  today's  production  rotorcraft. 
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Figure  21  -  Hub  lateral  vibration  comparison 

F.  GROUND  RESONANCE  ANALYSIS 
1 .  Frequency  Coalescence 

In  the  use  of  nonlinear  stiffness  to  avoid  ground  and 
air  resonance,  fundamental  questions  arise:  1)  In  what  way 
does  the  stiffness  effect  the  coalescence  of  the  regressing 


Base  Case  w/  Duffing  &  Incr  Omega 
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lag  mode  to  the  rigid  body  mode?  and  2)  Are  the  restoring 
moments  produced  of  reasonable  magnitude? 

Figure  22  is  the  Coleman  Plot  of  the  baseline  rotor 
with  and  without  Duffing  stiffness  included  at  the  blade 
root  end.  It  can  be  seen  from  the  plot  that  the  linear 
spring  only  case  has  coalescence  with  the  hub  lateral  mode 
at  a  rotor  operating  RPM  of  239,  as  depicted  earlier  in  the 
divergent  time  histories.  Figure  14  and  Figure  15.  With  the 
inclusion  of  the  stiffening  Duffing  spring,  the  non-rotating 
lag  natural  frequency  is  increased  due  to  the  increased  lag 
displacement.  This  blade  stiffening  shifts  the  regressing 
lag  mode  away  from  the  hub  modal  frequency  at  any  operating 
RPM  where  resonance  may  occur.  The  action  of  the  hardening 
Duffing  spring  keeps  the  frequency  coelesence  above  the 
instantaneous  rotor  RPM. 
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Fixed  System  Fan  Plot  With  Duffing  Stiffness 
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Figure  22  -  Coleman  plot  of  baseline  rotor  with  Duffing  stiffness 


It  should  be  noted  that  Duffing  springs  of  sufficiently 
low  magnitude  do  not  stabilize  the  rotor,  in  fact,  the  rotor 
RPM  that  results  in  ground  resonance  is  increased  slightly 
over  the  rotor  RPM  where  the  linear  rotor  would  experience 
ground  resonance.  For  this  reason,  the  Duffing  replacement 
of  auxilary  lag  damping  produces  less  hub  motion  for  rotors 
with  a  lag  frequency  in  the  upper  range,  0.6  or  0.7  per  rev, 
vice  the  lower  0.3  per  rev  range,  which  produces 
unacceptable  hub  and  blade  motion. 
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2.  Comparison  with  Theory  -  Baseline  Rotor 

A  Coleman  stability  plot  for  the  baseline  rotor 

follows.  This  plot  was  produced  from  Matlab®  code  that  was 
modified  from  Rafanello  [Ref.  9] .  One  can  see  the  region  of 
instability  between  225  RPM  and  266  RPM.  This  plot  is  for 
the  baseline  rotor  with  no  Duffing  stiffness  employed. 


Coleman  Stability  Plot  -  Baseline  Rotor 


Figure  23  -  Coleman  plot  for  baseline  rotor 

One  can  see  the  duplication  of  this  region  of 
instability  in  the  simulation  results.  Figure  24.  This 
stability  plot  matches  the  plot  from  Coleman  quite  well.  The 
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region  of  instability  begins  at  225  RPM  in  similar  fashion 
to  the  theory.  The  NPS  simulation  becomes  stable  at  261  RPM 
where  the  theory  does  not  predict  stability  until  266  RPM. 
This  difference  matches  the  results  found  by  both  Robinson 
and  Rafanello  in  that  Coleman's  linear  analysis,  as  applied 
by  Deutsch,  produces  slightly  conservative  stability 
criteria  for  ground  resonance  [Refs.  1,  9]. 


210  220  230  240  250  260  270 

Omega  (RPM) 

Figure  24  -  Baseline  rotor  region  of  instability 
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3 .  Comparison  with  Theory  -  Duffing  Rotor 

The  following  Coleman  stability  plot  is  for  the 
baseline  rotor  with  8  degrees  lag  amplitude  with  a  Duffing 
coefficient  20  times  the  linear  stiffness.  One  might  note 
that  Coleman's  is  a  linear  analysis.  This  linear  analysis  is 
applied  to  the  nonlinear  Duffing  rotor  at  a  single  instance 
in  time  with  a  specific  lead  angle,  and  is  not  valid  at  any 
other  blade  lead  angle.  One  Can  immediately  see  that  the 
nominal  rotor  speed  of  239  RPM  is  now  stable.  The  Duffing 
stiffness  in  the  rotor  continuously  drives  the  unstable 
region  above  the  RPM  of  the  rotor. 
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Coleman  Stability  Plot  -  Duffing  Rotor 
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Figure  25  -  Coleman  plot  for  Duffing  rotor 

4.  The  Region  o£  Instability  -  Duffing  Rotor 

The  overall  effects  of  including  Duffing  stiffness  in 
the  lead/ lag  dynamics  can  be  more  easily  discerned  in  Figure 
26.  Hub  lateral  stability  at  a  variety  of  rotor  speeds  is 
plotted  for  four  different  Duffing  blade  stiffnesses.  The 
introduction  of  Duffing  stiffness  at  10  times  the  linear 
stiffness  greatly  increases  overall  hub  motion  stability, 
but  does  not  stabilize  the  rotor.  A  Duffing  stiffness  just 
over  30  times  the  linear  stiffness  would  provide  a  stable 
baseline  rotor  for  all  rotor  speeds . 
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Figure  26  -  Change  in  region  of  instability  with  Duffing  stiffness  included 

It  should  be  noted  that  in  Figure  26  that  the  damping 
values  plotted  on  the  Y-axis  are  linear  approximations  for 
damping  obtained  from  a  nonlinear  system  time  history  (for 
all  but  the  Kd/K=0.0  line)  .  The  Hilbert  transform  method  was 
used  to  determine  damping  from  the  hub  lateral  time  history 
for  each  point  on  the  lines  depicted.  The  last  8  seconds  of 
a  9  second  simulation  were  used  as  inputs  to  the  damping 
analysis.  Example  time  histories  for  the  Kd/K=20.0  rotor 
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follow  in  Figure  27.  These  plots  were  taken  at  the  three 
separate  rotor  speeds  labeled. 


5.  Duffing  Rotor  Effect  on  Structural  Moments 

Figure  28  depicts  the  blade  restoring  moment  due  to  the 
combination  of  linear  and  nonlinear  stiffness  as  a  function 
of  blade  lag  displacement.  For  the  baseline  rotor  operated 
at  the  center  of  instability,  the  lag  amplitude  was 
approximately  8  degrees .  This  condition  increases  the 
restoring  moment  less  than  50%  from  the  linear  spring  only 
case. 


0  5  10  15 

Lead/Lag  Angle,  (deg) 

Figure  28  -  Lead/lag  restoring  moments 
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Keep  in  mind,  however,  that  a  blade  experiencing  these 
increased  loads  due  to  root  end  stiffness  encounters  no  lag 
moments  from  the  now  removed  auxiliary  lag  damper.  Also,  the 
8-degree  lead  angles  seen  in  the  baseline  case  with  Duffing 
stiffness  (Figure  16)  was  for  a  "worst  case"  at  the  center 
of  the  region  of  instability.  The  center  of  instability  is 
an  unlikely  nominal  design  RPM  for  a  rotor,  even  one  with 
nonlinear  stiffness.  For  the  cases  shown  with  increased  £2, 

the  blade  lead  angles  were  on  the  order  of  a  degree,  which 
shows  virtually  no  additional  restoring  moment  on  this  plot. 
The  addition  of  Duffing  stiffness  at  reasonable  rotor 
speeds,  therefore,  is  not  producing  blade  bending  moments 
that  current  rotor  designs  do  not  already  handle. 

6.  Stiff -in-Plane?  -  Duffing  Rotor 

In  ground  resonance  analysis,  one  additional  question 
arises:  Since  a  stiffening  term  is  being  employed,  is  the 
rotor  becoming  stiff-in-plane?  A  stiff-in-plane  rotor  is  one 
that  has  a  lead/ lag  natural  frequency  greater  than  the 
rotating  frequency  of  the  rotor  head.  For  the  baseline  rotor 
with  a  rotary  spring  Duffing  stiffness  coefficient  20  times 
the  linear  stiffness  coefficient,  the  rotor  only  goes  s tiff¬ 
in-plane  when  the  blade  lag  amplitude  exceeds  6  degrees  at 
239  RPM.  For  the  baseline  rotor  with  Duffing  stiffness 
employed,  the  rotor  can  go  stiff-in-plane.  If  the  rotor  is 
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operated  off  the  center  of  instability,  however,  it  does  not 
go  stiff-in-plane. 


Blade  Lead  Angle,  (deg) 


Figure  29  -  Soft  in-plane  frequencies  at  239  RPM 

As  depicted  in  Figure  30,  in  this  case,  the  rotor  goes 
stiff-in-plane  when  the  blade  lag  amplitude  exceeds  12 
degrees  for  this  baseline  rotor  with  Duffing  stiffness 
employed  at  310  RPM. 

At  the  increased  rotor  operating  speed  of  310  RPM,  the 
limit  cycle  lag  amplitude  was  less  than  a  degree.  For  lead 
angles  of  this  small  amplitude,  the  rotor  sees  very  little 
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migration  towards  stiff -in-plane  from  its  initial  lag 
frequency.  Since  the  exact  same  baseline  rotor  is  modeled, 
note  that  the  lag  natural  frequency  is  changed  due  to  the 
increase  in  rotor  RPM.  What  was  a  0.6  per  rev  lag  frequency 
rotor  is  now  a  0 . 6*  (239/310)  Q  =  0.46  per  rev  lag  frequency 
rotor. 


Lead/Lag  Frequencies  at  Omega  =  310  RPM 
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Figure  30  -  Soft  in-plane  frequencies  at  310  RPM 
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6. 


CHAOS  ANALYSIS 


1.  Response  to  Similar  Initial  Conditions 

Duffing' s  equation  is  well  documented  in  the 
literature.  In  Ueda  [Ref.  19],  parameter  combinations  that 
produce  chaotic  response  from  Duffing' s  equation  are 
summarized.  Since  chaotic  response  in  the  blade  lag  motion 
is  undesired,  a  check  of  the  simulation  time  histories  must 
be  completed  to  rule  out  chaotic  response  for  the  Duffing 
rotor  presented  in  this  paper. 

One  characteristic  of  chaotic  systems  is  vastly 
differing  time  histories  resulting  from  nearly  identical 
initial  conditions.  Figure  31  is  a  Poincare  section  of  the 
lateral  response  for  the  baseline  hub  with  Duffing 
stiffness.  Three  initial  conditions  of  hub  displacement  are 
plotted.  The  hub  response  depicted  in  this  plot  is  not 
chaotic;  the  sampled  points  on  the  phase  plane  are  nearly 
identical  for  the  three  cases.  The  sample  frequency  for  this 
Poincare  section  is  the  regressing  lag  frequency,  (Q-colag)  . 
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2.5 


Poincare  Map  of  Hub  Motion,  (Omegaregr) 
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Figure  31  -  Poincare  section  of  hub  lateral  motion 


2.  Poincare  Phase  Plane  and  Phase  Space  Analyses 

In  addition,  chaotic  motion  would  not  follow  a 
prescribed  form  in  phase  space.  Figure  32  is  the  blade 
lead/lag  motion  plotted  in  3-D  phase  space  with  the  blade 
lead  angle  rate  on  the  vertical  axis,  the  lead  angle  on  the 
radial  axis,  and  the  phase  plane  rotated  at  the  lead/lag 
natural  frequency.  The  initial  conditions  of  zero  lead  angle 
and  rate  quickly  begin  to  track  a  predictable  path  in  the 
phase  space,  showing  no  tendency  towards  chaos. 
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Figure  32  -  Phase  space  of  blade  lead  motion 


H.  DROOPING  ROTOR  STABILITY 

1.  Nonlinear  stiffness  provides  rotor  stability  for  a 
wide  range  of  physical  parameters 

Rotor  design  engineers  must  insure  rotor  stability  for 

a  wide  range  of  operating  parameters.  The  fuselage  rigid 

body  natural  frequency  in  roll  is  a  strong  function  of  the 
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fuselage  inertial  properties.  Different  fuel  loads,  cargo 
and  passenger  configurations  effect  fuselage  rigid  body- 
natural  frequencies.  Landing  gear  properties,  such  as  strut 
cleanliness,  strut  oleo  servicing,  and  tire  inflation  can 
also  effect  these  frequencies.  The  engineer  must  work  to 
obtain  a  stable  rotor  under  all  these  conditions.  In 
addition  to  fuselage  characteristics  changing,  rotor  blade 
properties  can  change.  Helicopter  operators  can  change  trim 
tab  and  pitch  link  settings,  effecting  the  way  the  blade 
tracks.  They  can  also  change  blade  tip  weights  in  the 
dynamic  balancing  process,  which  has  a  strong  effect  on 
blade  inertias . 

Historically,  larger  lead/ lag  dampers  than  were 
necessary  were  employed  in  final  designs  in  an  attempt  to 
maintain  rotor  stability  over  all  possible  operating 
conditions.  The  dampers  were  used  as  an  engineering  "safety 
net."  In  Figure  33,  the  Duffing  stiffness  in  the  rotor  is 
used  as  a  similar  safety  net  for  rotor  speed  reduction. 
Reduced  rotor  speeds  could  be  experienced  as  a  result  of 
system  failure  or  during  maintenance  checks,  in  addition  to 
rotor  engagement  and  disengagement.  This  figure  compares  a 
rotor  with  Duffing  stiffness  to  a  rotor  with  linear 
stiffness  only.  As  rotor  RPM  is  reduced,  both  rotors  enter 
the  region  of  instability  and  experience  resonance.  The 
linear  rotor  quickly  reaches  hub  displacement  amplitudes 
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that  could  cause  catastrophic  failure.  The  Duffing  rotor 
encounters  smaller  hub  displacements,  avoiding  loss  of  the 
helicopter.  In  this  case,  the  Duffing  stiffness  is  used  as  a 
rotor  safety  feature,  avoiding  the  large  hub  displacements 
that  result  in  failure  as  the  rotor  speed  is  reduced. 


Decreasing  Omega,  w/&  w/o  Duffing 
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Figure  33  -  Reduced  rotor  speed  comparison 


58 


IV.  LINEARLY  LINKED  ROTOR  SYSTEMS 


A.  INTRODUCTION 

The  helicopter  in  ground  resonance  could  be  compared  to 
a  simple  harmonic  oscillator.  The  rotor  acts  as  the  exciting 
force,  and  the  rigid  body  motion  of  the  helicopter  is  the 
response.  If  one  desires  to  avoid  resonance,  one  could  take 
one  of  two  approaches  to  the  problem.  Frequencies  of  the 
system  could  be  modified,  which  would  avoid  coalescence. 
Either  the  exciting  force  frequency  (the  rotor  regressing 
mode)  or  the  oscillator  natural  frequency  (the  fuselage 
rigid  body  natural  frequencies)  could  be  changed  to  avoid 
driving  the  system  into  resonance. 

The  other  approach  that  could  be  used  is  the 
modification  of  the  resonant  mode  shape.  Changing  the  mode 
shape  naturally  changes  the  mode's  resonant  frequency.  By 
linking  the  blades  together,  the  mode  shapes  of  the  blades' 
lead/ lag  motion  are  changed.  Modifying  the  mode  shape,  in 
turn,  alters  the  blade  natural  frequency  in  lead/ lag  and  the 
regressing  mode  of  the  rotor.  Therefore,  by  linking  the 
blades,  the  potential  exists  to  avoid  ground/air  resonance 
by  detuning  the  rotor-body  dynamics,  thus  avoiding  ground 
and  air  resonance  in  a  similar  manner  to  the  nonlinear 
stiffness  approach  already  discussed. 
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B.  ROTOR  WITH  LINKED  BLADES 

A  four  bladed  rotor  simulation  was  developed  similarly 
to  the  Coleman-like  three  bladed  rotor  used  in  the 
damperless  rotor  investigations  and  the  Comanche  five  bladed 
rotor  used  in  the  second  chapter.  This  baseline  rotor  was 
intentionally  chosen  to  be  unstable.  The  equations  of  motion 
for  this  new  four  bladed  simulation  were  also  derived  using 

Lagrange's  method  in  a  Maple®  worksheet.  Unlike  the  other 
rotor  simulations,  this  new  four  bladed  rotor  had  the 
additional  capability  of  linking  the  blades,  either 
elastically  or  rigidly.  In  addition  to  linear  stiffness  as  a 
function  of  the  angle  made  between  the  blade  and  its  nominal 
position  with  respect  to  the  hub,  in  this  case,  rotary 
stiffness  terms  were  added  as  functions  of  the  lead/lag 
angles  made  between  adjacent  blades.  All  linked  stiffness 
and  damping  values  were  linear. 

If  traditional  linear  root-end  stiffness  were  described 
as  Mn=k^*^n  for  the  nth  blade,  then  this  new  interlink 
stiffness  could  be  described  as: 


^inter(n.l)  (  Cn*l“  Cn  )  "’"^inter  (n)  ^Cn-1_ Cn  ^ 


Where  Kincer  is  the  interlink  stiffness  (moment/angle)  ,  and 
(Cn*i-C  J  and  (£n-i-C„)  are  the  instantaneous  lag  angle 
differences  between  the  nth  blade  and  its  two  adjacent 
neighboring  blades,  (n+1)  and  (n-1) .  The  simulation  was  also 
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given  the  capability  of  linking  the  blades  with  viscous 
dampers  which  were  functions  of  blade  lead  angle  rate 
differences : 

•M-n  =  Cinter(n+1)  (£n+l  —  Cn  )  +  Cinter(n)  (C  n-\  ~  £ n) 

Investigations  in  the  potential  payoffs  of  inter-blade 
linking  follow  in  the  next  three  sections .  The  baseline  case 
properties  for  the  four  bladed  interlink  rotor  may  be  found 
in  the  following  table: 

Table  3 


Q  =  25.0  rad/s 
(239  RPM) 

=  160.0  si 
(5150  lbs) 

“L.  =  4.0  si 
(129  lbs) 

e  =  0.5  ft 

R  =  10.5  ft 

®lag  /  &  =  0.6 

Linear  Damping  in  Lag  =  1.0% 

Linear  damping  in  Hub  Lateral  Motion  =2% 

Time  histories  for  the  baseline  case  blade  lead  angles 
may  be  found  in  Figure  34.  This  case  also  has  been  purposely 
set  inside  the  region  of  instability,  and  clearly  diverges. 
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Figure  34  •  Four  bladed  interlink  baseline 


C.  LINKED  BLADE  PARAMETRIC  STUDIES: 


1.  Linked  blade  &  hub  response  varying  blade 
interlink  stiffness  for  baseline  rotor 


The  effects  of  varying  blade  interlink  stiffness  for 


four  different  hub  damping  values  can  be  seen  in  Figure  35. 
For  the  cases  presented,  each  blade  is  elastically  linked  to 


each  of  the  two  neighboring  blades  in  series .  For  very  low 


interlink  stiffness,  KlinV/Kb1,^.  is  low,  the  simulations  are 
unstable,  similar  to  the  baseline  case,  as  would  be 


expected.  As  the  interlink  stiffness  is  increased,  however. 


Interlink  Baseline  Case,  239  RPM 
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a  sharp  increase  in  rotor  stability  occurs  in  the  region 
where  the  blade  interlink  stiffness  is  of  the  same  order  of 
magnitude  as  the  flexbeam  chord-wise  stiffness,  Kblade .  For 
these  cases,  previously  unstable  rotors  were  stabilized  for 
all  but  the  lowest  hub  damping  of  1%. 


Blade  Response  Damping  vs  Blade  Link  Stiffness  Ratio 


10"4  10'3  10‘2  10'1  10° 


K  /  K 
"link  '  "blade 

Figure  35  -  Blade  elastic  interlink  results 
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2.  Linked  blade  &  hub  response  varying  blade 
interlink  damping  for  baseline  rotor 

Interlinking  dampers  was  also  attempted.  All  four 
blades  were  interconnected  by  similar  dampers,  thus  putting 
the  interlink  dampers  in  series.  Recall  for  the  baseline 
four  bladed  rotor  that  linear  blade  damping  was  only  1.0%. 
For  inter-blade  damping,  the  same  four  hub  values  are 
plotted  in  Figure  36.  The  interlinking  of  blade  dampers  also 
proves  effective.  Unlike  the  case  of  the  interlinked 
stiffness,  this  increase  in  stability  is  to  be  expected 
since  the  overall  system  damping  has  been  increased,  even 
though  the  additional  dampers  are  implemented  in  less  than 
traditional  means. 
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Figure  36  -  Blade  damper  interlink  results 

3.  Linked  blade  &  hub  response  with  uneven  interlink 
stiffness  for  baseline  rotor 

Looking  at  unsymmetrical  blade  interlinking  was 

accomplished  so  as  to  change  the  relative  phasing  of  the 

blades  in  ground  resonance.  Note  that  for  the  divergent 

baseline  case  (no  linking) ,  each  blade  moves  90  degrees  out 

of  phase  with  its  neighbor  in  the  case  of  a  four-bladed 

rotor.  In  general  for  an  n-bladed  rotor,  Coleman  showed  that 


65 


during  divergence  the  blades  moved  at  a  phase  angle  of  2n/n 

(radians)  or  360/n  (degrees)  with  respect  to  each  other.  It 
was  expected  that  disrupting  this  phasing  would  change  the 
magnitude  of  the  rotor  CG  offset  from  the  center  of 
rotation,  in  turn  reducing  the  coupling  into  the  fuselage 
degrees  of  freedom.  This  is  due  to  the  fact  that  we  can 
change  the  natural  frequency  of  a  system  by  altering  its 
eigenvalue  (changing  the  spring  rate)  or  by  altering  its 
eigenvalue  (changing  the  mode  shape) .  In  Figure  37,  the  #1 
and  #2  blades  alone  are  linked.  This  action  does  reduce  the 
amplitude  of  the  blade  lead  angles  over  the  baseline 
unstable  case,  but  all  cases  using  this  method  diverged, 
nonetheless . 
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Blades  One  and  Two  Linked 


5  5.5  6  6.5  7  7.5  8  8.5  9 


Time  (sec) 

Figure  37  -  Blades  #1  and  #2  linked 

Linking  the  blades  in  pairs,  however,  showed  far 
different  results.  In  Figure  38,  the  #1  and  #4  blades  are 
linked,  and  the  #2  and  #3  blades  are  linked.  This 
interlinking  in  pairs  stabilizes  the  rotor.  For  a  linearly 
unstable  rotor,  blade  limit  cycle  amplitude  is  only  2 
degrees . 
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Blades  Four/One  and  Two/Three  Linked 


Time  (sec) 


Figure  38  -  Blades  #4&1  and  #2&3  linked 


By  linking  the  blades  in  adjacent  pairs,  the  four 
bladed  rotor  acts  more  like  a  two  bladed  rotor,  unable  to 
enter  ground  resonance.  This  method  of  stabilizing  the  rotor 
shows  promise  in  that  the  aerodynamic  and  handling  response 
characteristics  of  a  multi-blade  rotor  can  be  maintained, 
but  the  rotor  exhibits  the  mechanical  stability  normally 
found  in  two  bladed  rotors . 
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4.  Interlink  Dampers  Versus  Conventional  Dampers 

In  Figure  39,  traditional  blade  dampers  are  compared  to 
interlink  dampers.  Essentially,  the  four  lead/lag  dampers 
that  might  be  found  on  a  rotor  have  been  detached  from  the 
hub  and  remounted  between  the  four  blade  roots .  Identical 
damping  values  were  used  for  the  two  lines  plotted. 


Blade  Response  vs  Auxiliary  Damper  Strength 


1 1 _ i - 1 _ i _ i _ i _ i _ i _ i _ I 

'1  23456789  10 


Aux.  Damper  Damping  (%) 

Figure  39  -  Conventional  dampers  versus  interlink  dampers 

The  interlink  damping  proves  more  effective  for  a  given 
damping  value.  Instead  of  the  dampers  acting  independently 
for  each  blade,  the  new  system  takes  advantage  of  the  fact 
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that  the  blades  are  90  degrees  out  of  phase  with  each  other, 
essentially  amplifying  the  input  to  the  dampers.  Since 
weight  is  always  a  driving  issue  in  aircraft  design,  if  the 
dampers  were  mounted  between  the  blades,  instead  of  from  the 
blade  to  the  hub,  smaller  dampers  could  be  installed  and 
still  maintain  the  desired  damping  in  the  blade  chord-wise 
dynamics . 
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V.  ROTOR  WITH  UNEVEN  BLADE  SPACING 

As  long  as  unsyxnmetrical  rotors  are  being  investigated, 
the  effect  of  modifying  blade  spacing  was  also  attempted. 
Blade  motion  time  histories  from  a  60-120-60-120  degree 
spread  rotor  can  be  found  in  Figure  40.  This  rotor  layout  is 
similar  to  the  AH-64  Apache  tail  rotor,  where  the  four 
blades  are  not  positioned  90  degrees  apart,  but  instead 
arranged  in  an  X  pattern,  not  a  perfect  cross.  The  Apache¬ 
like  rotor  also  exhibited  blade  divergence.  The  blade  pairs 
diverged  at  different  rates,  but  the  same  overall  response 
was  noted:  mechanical  instability  was  not  avoided  for  any  of 
the  rotor  speeds  investigated  for  the  case  where  uneven 
blade  spacing  was  introduced. 
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Baseline  Rotor  With  60-120-60-120  Blade  Spread 


123456789 
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Figure  40  -  Baseline  rotor  with  60-120-60-120  blade  spread 


72 


VI.  GROUND  AND  AIR  RESONANCE  ANIMATION 

Obviously,  numerous  cases  and  configurations  were 
simulated  in  the  course  of  this  research.  All  these 
simulations  were  performed  on  a  PC  running  Simulink®. 

Utilizing  the  numerical  analysis  capabilities  of  Matlab® 

including  the  Hilbert  damping  analysis,  numerous 
quantitative  parameters  were  generated. 

In  an  attempt  to  provide  a  better  qualitative  picture 
for  what  the  rotor  was  doing,  an  animation  routine  was 
written  that  would  show  the  engineer  visually  what  the 
blades  and  hub  were  doing  in  the  horizontal  plane.  The 
animation  routine  written  as  part  of  this  research  plots  the 
rotor  blades  as  simple  lines  at  the  computed  lead  angle  from 
the  simulation  at  every  time  step.  The  depicted  rotor  does 
not  rotate,  so  that  blade  action  can  be  observed.  The  center 
of  rotation  is  also  changed.  The  animated  blade  center  is 
plotted  at  the  hub  X  and  Y  position  with  respect  to  the 
inertial  frame,  also  obtained  from  the  simulation  at  every 
time  step.  In  this  way,  the  outward  spiral  of  the  hub  is 
also  depicted  as  part  of  the  animation. 

Though  no  numerical  results  are  obtained,  the  animation 
routine  was  found  very  useful  in  providing  an  overview  of 
rotor  states.  Four  example  frames  of  this  simulation  are 
shown  in  Figure  41.  These  frames  show  the  progression  of 
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blade  movement  in  a  typical  ground  resonance  scenario. 
Notice  how  the  animation  illustrates  the  spiral  motion  of 
the  CG  of  the  rotor  around  the  center  of  rotation.  The  90 
degree  blade  phasing  is  also  easily  discernible. 
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Rotor  Simulation  Animation 

Simulation  time: 

9.9167  sec 

©  *  t  =  60°  in  Blade  #1  Lead  Cycle 
(blade  *1  is  at  6  o’clock  in  figure) 


9.9833  sec 

©  *  t  =  120°  in  Blade  *1  Lead  Cycle 


10.1167  sec 

©  *  t  =  240°  in  Blade  #1  Lead  Cycle 


10.1833  sec 

©  *  t  =  300°  in  Blade  ul  Lead  Cycle 


Figure  41  -  Rotor  animation 
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VII.  DAMPERLESS  ROTOR  DEVELOPMENT  CONCLUSIONS 

With  the  advent  of  hingeless  and  bearingless  main  rotor 
helicopters,  interest  in  nonlinear  rotor  dynamics  has 
increased  significantly  in  recent  years  as  an  area  to 
research  for  improved  rotor  stability.  With  the  rapid 
evolution  of  high-speed  processors  and  recent  developments 
in  software  technology,  we  now  have  the  tools  available  to 
explore  rotor  system  nonlinear  dynamics  and  potentially 
achieve  helicopter  rotor  mechanical  stability. 

This  dissertation  furthered  the  development  of  a  non¬ 
linear  analysis  tool  that  has  the  capability  of  modeling  a 
rotor's  nonlinear  material  properties.  Most  importantly,  the 
simulation  did  not  incorporate  assumptions  or  ordering 
schemes  that  limit  the  number  of  terms  in  the  equations  of 
motion  and  subsequently  skew  the  output.  In  the  simulation 
used,  not  a  single  term  was  eliminated. 

Potential  methods  for  passively  controlling  lead/lag 
motion  of  helicopter  blades  have  been  developed.  One  of  the 
methods  presented  utilizes  Duffing-type  nonlinear  stiffness 
in  the  root  end  of  the  rotor  blade  in  addition  to  linear 
spring  and  damping  properties.  The  Duffing  spring  provides  a 
restoring  force  that  is  a  cubic  function  of  the  blade 
lead/lag  deflection.  This  nonlinearity  shifts  the  first  in¬ 
plane  natural  frequency  of  the  blade  as  a  function  of  the 
blade  lag  angle.  This  shift  in  blade  natural  frequency 
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effectively  varies  the  regressing  mode  frequency.  The 
resulting  shift  in  the  regressing  mode  prevents  coalescence 
with  hub  resonant  frequencies.  For  a  linearly  unstable 
rotor-fuselage  system,  inclusion  of  the  Duffing  spring 
results  in  limit  cycle  behavior  in  the  blade  lead/ lag 
degrees  of  freedom  and  in  hub  translational  motion. 

An  analysis  of  the  limit  cycle  motion  effect  on  hub 
vibrations  was  also  completed.  The  replacement  of  blade 
auxiliary  lag  dampers  with  nonlinear  springs  results  in 
increased  hub  vibration  levels  if  the  resulting  system  is 
highly  unstable  in  the  linear  analysis .  Rotors  that  are 
marginally  stable,  or  stable,  in  the  linear  analysis  produce 
vibration  levels  that  are  not  uncommon  in  production 
rotorcraf t . 

The  primary  conclusions  of  this  research  on  damperless 
rotors  are: 

•  Application  of  nonlinear  stiffness  properties  to  the 

blade  lead/lag  degrees  of  freedom  can  effectively  improve 
helicopter  ground/air  resonance  stability 

characteristics . 

•  Stability  characteristics  are  improved  through  frequency 
mismatch,  not  by  dominating  system  dynamics  with  high 
damping . 

•  Blades  with  low  lead/lag  natural  frequencies  have  greater 
potential  for  producing  unacceptable  vibration  levels  in 
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both  dynamic  and  non-rotating  components  if  Duffing 
stiffness  is  employed. 

•  Interblade  rigid  and  elastic  coupling  has  the  potential 
to  eliminate  helicopter  unbounded  blade  motion. 

•  In  this  case,  the  interblade  coupling  achieves  stability 
by  altering  the  phase  of  relative  blade  response.  That 
is,  the  response  eigenvector  has  been  altered. 

•  Interblade  viscous  damping,  selectively  applied,  can  be  a 
more  efficient  use  of  blade  auxiliary  dampers  than  the 
conventional  blade-to-hub  arrangement. 

•  Uneven  blade  spacing,  for  the  case  examined  (60-120-60- 
120),  did  not  prove  effective  in  avoiding  helicopter 
ground  and  air  resonance. 
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APPENDIX  A.  COMANCHE  FROUDE  SCALE  MODEL  SIMULATION 


Comanche  Rotor  Design  Overview 

Engineering  data  on  the  Boeing/ Sikorsky  RAH- 6 6  Comanche 
Froude  scale  wind  tunnel  model  was  obtained  from  Sikorsky 
Aircraft  as  part  of  a  research  program  initiated  through  the 
National  Rotorcraft  Technology  Center  (NRTC)  [Refs.  Al,  A2] . 
Additional  material  on  the  Comanche  rotor  is  provided  by 
Tar z an in  and  Panda  [Ref.  A3] . 

The  Boeing/Sikorsky  report  obtained  by  NPS  outlines 
testing  procedures,  results,  and  basic  rotor  design  geometry 
of  the  Comanche  model  rotor.  The  Comanche  rotor  design 
utilizes  a  bearingless  main  rotor  (BMR)  flexbeam  arrangement 
for  blade  lead/lag,  flap,  and  feathering  degrees  of  freedom. 
A  cuff,  which  is  stiff  in  torsion  and  bending,  encases  the 
flexbeam  and  provides  input  torques  in  the  feathering  axis 
to  the  blade  at  the  outboard  end  of  the  flexbeam,  at  the 
cuff /blade  joint.  A  damping  device  called  the  snubber -damper 
is  mounted  between  the  inboard  end  of  the  cuff  and  the  hub 
to  provide  additional  damping  and  stiffness  in  the  blade 
lead/ lag  degrees  of  freedom.  A  depiction  of  this  rotor  is 
given  in  Figure  Al,  reproduced  from  Panda. 
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Figure  A1  -  Froude  scale  Comanche  rotor  design  [From  Ref.  A3] 


Description  of  Nonlinear  Snubber-Damper 

The  snubber-damper  is  made  up  of  alternating  layers  of 
metal  and  elastomer,  which  provides  damping  when  the  layers 
are  put  in  shear.  Unfortunately,  the  damping  and  stiffness 
values  for  the  snubber-damper  are  not  constant;  they  are  a 
strong  function  of  the  snubber-damper  displacement 
amplitude.  The  initial  snubber-damper  designs  provided  very 
little  damping  for  small  deflections. 
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In  addition  to  the  low  damping  values,  snubber-damper 
stiffness  tended  to  rise  sharply  for  very  small  amplitudes. 
Plots  of  the  snubber-damper  stiffness  and  loss  factor  were 
published  by  Panda  and  follow  as  Figure  A2,  [Ref.  A3]  .  The 
"Elastomeric  Damper"  depicted  is  the  device  used  in  the 
1992-1993  wind  tunnel  tests  and  the  "Fluidlastic®"  damper 
was  the  result  of  a  redesign,  tested  in  1995. 
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Figure  A2  -  Snubber-damper  stiffness  and  loss  factor  [From  Ref.  A3] 

Inherent  in  the  design  of  these  devices,  the  snubber- 
dampers  exhibit  hysteresis  in  their  blade  restoring  force  as 
a  function  of  displacement.  Panda  reported  the  Comanche 
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elastomeric  snubber-damper  as  having  sufficiently  low 
damping  at  small  displacements  that  this  device  caused  limit 
cycle  oscillations  in  the  blade  chord-wise  motion.  This 
limit  cycle  motion  produced  undesired  vibrations  in  the 

coupled  rotor- fuselage  system.  The  Fluidlastic®  snubber- 
damper  alleviated  the  severity  of  the  limit  cycle  motion 
problem.  One  can  see  less  reduction  in  the  Fluidlastic® 

damper  loss  factor  at  small  displacement  amplitudes  in  the 
lower  plot  of  the  previous  Figure. 

Flexbeam  and  Cuff  Modeling 

In  order  to  accurately  model  the  1/6  scale  Comanche 
wind  tunnel  test  rotor,  flexibility  in  the  blade  and 
flexbeam  had  to  be  accurately  modeled.  Initially,  a  linear 
beam  approximation  was  used.  This  was  done  primarily  as  a 
first  iteration  in  the  model  design  process.  The  linear 
theory  also  provided  data  for  verification  of  a  Myklestad 
method  that  was  coded  in  Matlab®  as  an  additional  task  in 
this  research. 
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Comanche  Model  Assembly 

The  Comanche  1/6  scale  flexbeam  and  blade  were  modeled 
using  the  Myklestad  program.  Example  output  of  a  1-inch 
blade  tip  deflection  in  lag  follows  as  Figure  A3 . 


Myklestad  Method 


Radial  position  (inches) 


Figure  A3  -  Combined  flexbeam/blade  analysis 

Using  the  Myklestad  program,  the  equivalent  offset  was 
found  to  be  a  constant  value  of  8.52  inches  from  the  center 
of  rotation,  with  the  snubber-damper  not  installed.  This 
equivalent  offset  is  highlighted  in  Figure  A4  by  a  small 
circle . 
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Equivalent  Offset  by  Myklestad  Method 


0  5  10  15  20  25  30  35  40 


Radial  location  (in) 

Figure  A4  -  Flexbeam/blade  deflection 

From  Figure  A4,  it  can  be  easily  seen  that  the  outboard 
aerodynamic  portion  of  the  blade  moves  as  a  rigid  body  for 
in-plane  bending  at  the  frequencies  of  interest.  Therefore, 
the  portion  of  the  blade  outboard  of  the  equivalent  offset 
was  modeled  as  a  single  degree  of  freedom.  Since  the  slope 
change  in  the  combined  model  is  at  the  root  end  in  the  soft 
region  of  the  flexbeam,  an  articulated  type  attachment  could 
be  used  for  mathematical  modeling  of  the  hingeless  rotor  in 
the  Simulink®  simulation. 
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Using  the  blade  mass  distribution  data  provided,  a 
blade  mass  and  moment  arm  were  computed.  A  blade  moment  of 
inertia  was  calculated  from  this  data  about  the  equivalent 
offset.  This  inertia  was  then  broken  down  into  an  equivalent 
mass,  ir^,  and  moment  arm,  1.  The  blade  was  modeled  as  a 
point  mass  at  a  distance  1  from  the  offset,  such  that  the 
total  mass  of  tho  actual  blade  matched  the  point  mass  value 
in  the  model.  For  accurate  modeling  of  the  lead/lag 
dynamics,  it  was  desired  that  e*S/I  for  the  actual  blade 
match  e/1  for  the  point  mass  model,  where  S  is  the  first 
mass  moment  of  the  blade  and  I  equals  the  second  mass  moment 
of  the  blade.  The  blade  model  was  determined  to  be  a  0.0129 
slug  point  mass  at  a  radial  location  of  30.06  inches. 

Model  parameters  for  the  portion  of  the  blade  outboard 
of  the  effective  offset  are  compared  to  the  actual  blade  in 
the  following  table: 


Table  A1 


Froude  Model 

NPS  Model 

Blade/cuff  moment  of  inertia,  I 

4.5212  sl*inA2 

Blade/cuff  mass  moment,  S 

0.2099  sl*in 

0.0129  si 

I/S 

21.5407  in 

Blade  lag  offset 

8.5435  in 

Point  mass  radial  location 

30.0642  in 

The  cuff  is  far  stiffer  than  the  flexbeam,  providing  an 
essentially  rigid  sleeve  around  the  flexbeam.  For  lead/ lag 
motion,  the  flexbeam  simply  bends  inside  the  interior 
dimensions  of  the  cuff.  Since  the  snubber -damper  attaches  at 
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the  inboard  end  of  the  cuff,  the  snubber-damper 
displacements  were  considered  to  be  a  function  of  an 
extension  of  the  cuff /blade  joint  angular  displacement  and 
off-axis  deflection,  as  seen  below  in  Figure  A5: 


Figure  A5  -  Snubber-damper  motion  as  a  function  of  cuff/blade  joint  motion 

Assuming  a  constant  equivalent  offset  for  blade 
lead/ lag  motion,  a  geometric  gain  was  computed  for  snubber- 
damper  displacement  as  a  function  of  blade  lead/lag  angle. 
Using  the  physical  properties  of  the  blade  root  end,  the 
geometric  gain  was  computed  to  be  5.937  inches  of  snubber- 
damper  deflection  in  shear  for  every  radian  of  blade  lead, 
assuming  small  lead/lag  angles. 

Snubber -Damper  Model 
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The  NPS  nonlinear  rotor  simulation  was  modified  to 
incorporate  nonlinear  snubber-damper  coefficients  in  both 
stiffness  and  damping.  The  loss  factors  made  available  in 
the  Boeing/Sikorsky  report  were  first  converted  to 
dimensional  damping  coefficients.  The  loss  factor  is  defined 
as  the  out-of-phase  force  divided  by  the  in-phase  force. 
Since  the  experimental  loss  factor  data  in  the  report  was 
recorded  at  a  fixed  driving  frequency  of  10  Hz,  the  loss 
factor  data  could  be  converted  to  dimensional  damping 
values.  The  dimensional  damping  value  may  be  expressed  as: 


c(x)  = 


yk{x) 

0) 


where  T|  is  the  loss  factor,  k  is  the  nonlinear  snubber- 


damper  stiffness,  CO  is  the  forcing  function  frequency,  and 
x  is  the  snubber-damper  displacement. 

Traditional  hydraulic  and  viscous  damping  methods  were 
attempted  in  the  nonlinear  snubber-damper  mathematical 
model,  in  addition  to  higher  order  curve  fits.  A  fourth 
order  fit,  c  (x)  =cix‘+c3x +c2x2+c1x+c0,  in  displacement  was 
decided  upon  as  the  best  trade-off  between  accuracy  and 
model  order  reduction.  A  linear  relationship  was  used  for 
snubber-damper  displacements  greater  than  0.6  inches. 

The  polynomials  formulated  for  both  the  snubber -damper 
damping  and  stiffness  formulations  follow: 
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Table  A2 


1992  Elastomeric  snubber-damper: 

c=-842875x‘+87576 . 8x -2847 . 07xz+25 . 01x+0 . 47 
k=3  3  62410  84x*-3  8  666  827xJ+1547  897x^-26483 . 84x+233 .32 
1995  Fluidlastic®  snubber-damper; 

c=116695 . 68x'-14022 . 78xJ+609  .  Ox"- 12 . 05x+0 . 59 
k=38512436 . 84x*-4491019 . 58xJ+18333 6 . 70x^-3147 . 25x+79 .71 


where  c  is  in  lbs/in/s,  k  in  lbs/in,  and  x  in  inches.  The 
stiffness  and  damping  results  are  compared  to  the  snubber- 
damper  laboratory  data  provided  in  the  plots  that  follow: 
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Damping  (ib/ips) 


Figure  A6  -  Snubber-damper  stiffness  and  damping  properties 
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NPS  Simulation  Results  for  the  1/6  Scale  Comanche 

Example  results  of  the  NPS  rotor  simulation  follow.  The 
inputs  for  these  plots  may  be  found  in  the  following  table: 


Table  A3  -  Helicopter  physical  and  aerodynamic  parameters 


Helo  Physical  and  Aerodynamic  Parameters  | 

%  Rotor  speed,  rad/sec 

0mega=870* (pi/30) 

%  Equivalent  lag  offset,  inches 

el=8 . 5235 

%  Equivalent  blade  length,  inches 

1=21.5407 

%  Snubber  radial  location,  inches 

Rs=2 . 5984 

%  Lead/lag  stop  position,  rads 

No  lag  stops  modeled 

%  Mass  of  rotor  blades,  slugs 

mb=0 . 0129 

%  Effective  mass  of  fuselage,  slugs 

Mx=  0.6452 

My=  1.3802 

%  Blade  azimuth  phase  angles,  rads 

Phi ( 1 : 5 ) = [ 0  246  8]*pi/5 

%  geometric  gain  of  snubber 
movement  from  lag,  in/rad 

gg=5  o 1121 

%  Spring  stiffness  polynomial  for 

Kbeam= [0,0, 0,0 ,660. 9141] 

lead-lag,  in* lbs /radian)  (1995 

Ksnub= [38512436 . 84*ggA4, 

snubber) 

-4491019. 58*ggA3, 18333 6. 70*ggA2, 
-3147 . 25*gg, 79.71] 

* (el-Rs) A2 

Kpoly=Kbeairw-Ksnub 

%  Lead-lag  stop  spring 

Ks=0 

constants,  in* lbs /radian 

No  lag  stops  modeled 

%  Damping  polynomial  in  lead- lag 

Cpoly= [116695 . 68*gg^4 , 

(in*lb/ (rad/sec) )  (1995  snubber) 

-14022 . 78*gg^3 , 609 . 0*gg^2 , 

-12 . 05*gg, 0.59]  *(el-Rs)A2 

Linear  springs  in  translation. 

K(l) =297 . 38 

lbs /in 

K (2 ) =217 . 93 

Linear  damping  in  translation. 

c (1) =1.29 

(lbs/ (in/ sec) ) 

c (2 ) =1 . 88 

Non-linear  damping  in  translation, 
lbs/ (in/sec) A2) 

v=0 

Fuselage  initial  displacements,  in 

xXi=0.1 

(all  other  IC's  =  0.0) 

xYi=0 . 1 

The  fourth  order  snubber-damper  model  is  employed  in 
the  simulation  depicted  in  the  following  time  histories.  All 
blades  are  given  initial  conditions  of  zero  lead  angle 
displacement  and  zero  velocity.  All  rotor  and  hub  dynamics 
are  initiated  by  a  0.1  foot  displacement  in  both  the  x  and  y 
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dimensions  of  hub  motion.  Hub  response  from  these  initial 
conditions  is  shown  in  Figure  A7.  This  initial  condition  is 
equivalent  to  a  lateral  restraining  force  on  the  hub  of 
36.87  lbs.  being  released  at  time  t=0.0  in  the  simulations. 
This  comparatively  large  impulse  was  set  to  insure  coupling 
of  the  rotor /body  dynamics . 


Hub  Translational  Time  Histories 


Figure  A7  -  Hub  translational  time  histories 

The  Comanche  wind  tunnel  tests  utilized  a  shaker  in  the 
fuselage  to  initiate  rotor/body  coupling  prior  to  recording 
time  histories  for  analysis.  Since  the  initial  conditions 
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varied  for  each  test  in  the  tunnel,  duplication  of  the  exact 
initial  conditions  was  impossible  with  the  simulation,  and 
the  fixed  initial  conditions  on  the  hub  displacement  were 
used.  The  Boeing /Sikorsky  report  describes  this  process  as 
follows : 

"The  blade  lead-lag  mode  was  excited  by  the 
swashplate  cyclic  shaker  operating  at  the  fixed 
system  lead-lag  regressing  frequency.  When 
required,  a  fine  adjustment  was  made  to  maximize 
the  lead-lag  response  of  the  model.  The  amplitude 
was  monitored  (on  line)  by  a  spectrum  analyzer. 
When  the  forced  response  was  judged  to  be  maximum, 
the  shaking  was  stopped,  marking  the  beginning  of 
the  transient  decay.  The  following  10  seconds  of 
chord  bending  and  fuselage  roll-pitch  data  was 
collected  and  analyzed  by  a  moving-block  method  to 
determine  the  frequency  and  damping  ratio  of  the 
mode  being  investigated."  [Ref.  A2] 


Response  to  the  relativey  stiff  rotor  blades  is  shown 
in  Figure  A8 .  The  Comanche  lead/ lag  frequency  was  reported 
to  be  0.71Q0  for  the  1995  system,  0.692 Q0  for  the  1992 
system. 
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Blade  Lead/Lag  Motion  Time  Histories 


Time  (s) 


Figure  A8  -  Blade  lead  angle  time  histories 

Damping  for  each  of  the  five  blades  was  determined 
using  a  Hilbert  transform  method  [Ref.  A5]  .  The  damping 
value  for  the  blades  was  then  averaged  for  the  final  8 
seconds  of  the  9  second  time  history.  This  average  was  then 
converted  to  a  fixed  system  damping  value  replicating  the 
method  used  in  the  Comanche  test: 

"The  rotating  system  damping  ratio  is  converted  to 
a  fixed-system  damping  ratio  by  multiplying  by  the 
ratio  of  the  rotating-system  lead-lag  mode  frequency 
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(©£)  to  the  fixed-system  lead-lag  regressing  mode 
frequency  (Q-co^)  . "  [Ref.  A2] 

An  example  output  of  the  Hilbert  damping  analysis  for 
the  #5  blade  follows: 


Hilbert  Method  Damping  Determination 


123456789 


Time  (sec) 


0  100  200  300  400  500  600  700  800  900 


Freq,  rad/s 

Figure  A9  -  Blade  motion  damping  analysis 

The  top  subplot  is  the  time  history  of  blade  lead  angle 
that  was  analyzed.  The  second  subplot  is  the  natural 
logarthm  of  the  absloute  value  of  the  Hilbert  transform  of 
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the  signal.  A  first  order  fit  to  this  plot  is  also  shown. 
The  damping  of  the  signal  was  then  computed  by  dividing  the 
negative  of  the  slope  by  the  time  history  frequency.  The 
primary  frequency  of  the  time  history  is  determined  by  power 
spectral  density  using  Matlab®'s  discrete  fourier  transform 
function,  "fft.m."  The  results  of  the  power  spectral  density 
are  shown  in  the  third  subplot.  This  Hilbert  transform 
method  works  well  for  signals  dominated  by  a  single  mode. 
Smith  obtained  more  accurate  results  with  the  Hilbert 
transform  method  than  with  Moving  Block  for  signals  with  a 
single  mode  present  [Ref.  A5] .  Results  of  the  Hilbert 
transform  damping  analysis  of  the  two  hub  degrees  of  freedom 
follow  in  Figures  A10  and  All. 
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Hub  X  (in) 


Hilbert  Method  Damping  Determination 
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Figure  A10  -  Hub  X  direction  damping 


Hilbert  Method  Damping  Determination 
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Figure  All  -  Hub  Y  direction  damping 


The  frequencies  obtained  from  the  NPS  simulation  are 
summarized  in  the  following  table: 

Table  A4 


NPS  Simulation 

Boeing/ Sikorsky  report 

Hub  X 

3.23  Hz 

3.20  Hz 

Hub  Y 

1.99  Hz 

2.00  Hz 

Lead /Lag 

0.76Q 

0.692Q0 
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Comparison  with  Test  Data 

The  goal  of  the  application  of  the  NPS  simulation  to 
the  Comanche  scale  rotor  data  was  the  matching  of  fixed 
system  damping  values  to  the  test  data  and  to  UMARC.  The 
plot  that  follows  reproduces  the  Boeing/ Sikorsky  report's 
Figure  50b,  from  the  1992  test  results  [Ref.  Al]  .  The  1992 
test  data  was  chosen  as  a  more  rigorous  test  of  the  NPS 
simulation,  since  the  snubber -damper  used  in  the  1992  tests 
exhibited  greater  nonlinear  behavior.  The  first  plot.  Figure 
A12,  shows  the  NPS  simulation's  initial  results: 


Replication  of  Boeing  Sikorsky  Figure  50b,  1992 


700  750  800  850  900  950  1000 

Omega,  RPM 

Figure  A12  -  Initial  NPS  simulation  results 
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The  geometric  gain,  the  snubber -damper  shear  to  blade 
lead  angle  ratio,  was  considered  a  soft  parameter  and  varied 
in  an  atttempt  to  obtain  more  accurate  matching  with  the 
existing  data. 

As  the  geometric  gain  was  increased,  slightly  better 
results  were  obtained.  The  following  plot  shows  the  results 
of  increasing  the  geometric  gain  by  a  factor  of  1.5.  From 
these  results,  the  wind  tunnel  model,  it  is  assumed,  was 
seeing  more  snubber-damper  movement  than  the  modeling  method 
assumptions  had  called  for.  This  soft  parameter,  the 
geometric  gain,  is  related  to  the  multiple  load  path  problem 
alluded  to  earlier.  In  brief,  the  basic  assumptions 
underestimated  the  amount  of  snubber-damper  motion  with  lead 
angle  changes.  Recall  only  lead/lag  motion  is  modeled,  and 
the  bearingless  Comanch  rotor  model  experienced  coupling  in 
the  lead/lag  dynamics  from  other  model  degrees  of  freedom. 

With  the  advantage  of  having  the  test  data  results  a 
priori,  the  model  was  tuned  to  best  replicate  the  wind 
tunnel  data  using  the  geometric  gain  as  a  tunable  parameter. 
With  the  higher  geometric  gain  mentioned,  the  NPS  simulation 
appeared  to  replicate  the  model  data  slightly  better  than 
UMARC  at  the  lower  RPM's  tested.  At  higher  rotor  speeds,  the 
difference  between  the  two  simulations  is  negligible. 
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Fixed  System  Damping  (%critical) 


Replication  of  Boeing  Sikorsky  Figure  50b,  1992 


Omega,  RPM 


Figure  A13  -  Final  NPS  simulation  results  with  increased  geometric  gain 


Comanche  Modeling  Conclusions 


•  If  the  physical  rotor 

exhibits 

migration  of 

the 

effective 

offset  with 

varying 

lead  angles 

and 

velocities. 

forcing  a 

hingeless 

rotor  into 

an 

articulated 

model  contaminates  the 

results . 

•  The  inclusion  of  aerodynamics,  coupling  in  the  rotor 
lead/ lag,  flap,  and  torsion  should  be  included  in 
the  modeling  of  a  bearingless  rotor. 
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•  Modeling  the  blade  portion  that  is  outboard  of  the 
effective  offset  with  only  two  parameters,  blade 
length  and  tip  mass,  is  incapable  of  matching  all 
the  inertial  properties  required  to  accurately 
simulate  both  lead/lag  and  hub  motion. 

•  Additional  degrees  of  freedom  are  required  in  the 
lead/lag  dynamics  to  better  approximate  multiple 
load  paths  from  the  blade  to  the  hub. 

•  For  a  fledgling  code,  the  NPS  simulation  shows 
potential  in  modeling  rotors  with  nonlinear  physical 
parameters  and  further  development  should  be 
considered. 
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APPENDIX  B.  MYKLESTAD  BEAM  MODELING  METHOD 


In  1944,  Myklestad  published  a  method  for  finite 
element  analysis  of  beams  in  bending  for  vibrations 
applications  in  his  book,  "Vibration  Analysis"  [Ref  Bl] .  His 
method  breaks  a  beam  into  segments  that  have  essentially 
linear  behavior  over  the  segment  length.  Equations  of  motion 
of  each  of  the  separate  elements  are  then  tied  together  as 
separate  degrees  of  freedom.  In  Myklestad' s  method,  each 
beam  element  is  defined  as  a  point  mass  and  a  massless 
elastic  element  in  bending,  as  depicted  in  Figure  Bl, 
modified  from  Gerstenberger  and  Wood  [Ref.  B2] . 


Figure  Bl  -  Typical  Myklestad  beam  element  [After  Ref.  B2] 
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By  defining  element  members  in  this  way,  a  beam  with 
inertial  and  mass  properties  that  change  as  a  function  of 
beam  longitudinal  position  may  be  well  approximated.  Shear, 
moment,  angular  displacement  and  beam  deflection  may  be 
described  for  each  element  as  follows: 


Sn+\=Sn+mnynQ)2 
^ n+ 1  ^ n  ^n,n+l  *^n+l 


0n+1  =<=>„ 

1  "  EIn  n+1  2 EnIn  n+I 

yM -y.  +4,.,©.  ~kfrs^ 

n  n  n 


where  to  is  the  first  beam  natural  frequency  in  bending,  lnntl 

is  the  segment  length,  En  is  the  element  stiffness,  and  In  is 
the  element  second  area  moment. 

Myklestad's  method  may  be  modified  from  the  simple  beam 
in  bending  to  rotating  beams  by  the  inclusion  of  centrifugal 
and  blade  tension  terms.  For  this  case,  •  the  x-y  plane  is 
transformed  to  a  radial  and  off-axis  displacement  frame,  as 
depicted  in  Figure  B2,  also  reproduced  from  Gerstenberger 
and  Wood. 
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Figure  B2  -  Typical  rotating  Myklestad  beam  element  [From  Ref.  B2] 

The  equations  for  a  Myklestad  beam  element  in  the 
rotating  frame  follow: 

Sn  =Sn+i+mn(6)2  +Q2)xn 

Tn  =Tn+l+mnQ2rn 

Mn  =Mn+l+SnJnn+1  ~Tn+l(xn+i  -X„) 

=  ®n+l  —  ^n+l^n,n+l )  +  ^  n+l^n,n+i  —  ^n+l^n,n+l 

Xn  =  -^n+l  +®Jn,n+ 1  +  ^n+l®n+l^n,n+l  —  ^ n+l^ n,n+l  ~  ^ n+\^n,n+] 

where  Vn  n+l  —lntn+\/EnIn,  U n  n+l  —ln<n+l  /2EnIn,  and  Gn  n+]  =  ln  n+\  /2>Enln . 

A  program  written  in  the  Matlab®  application  language 
that  analyzes  a  beam  using  Myklestad' s  method  may  be  found 
elsewhere  in  the  Appendicies.  This  program  analyzes  beam  or 
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rotor  blade  inertial,  geometric,  and  stiffness  data  for  any 
number  of  beam  segments,  up  to  the  limitations  of  the 
Matlab  host  program  used.  Both  rotating  and  non-rotating 
beams  can  be  analyzed. 

The  program  also  has  the  capability  of  applying 
different  boundary  conditions  at  the  ends  of  the  beam.  An 
articulated  blade  requires  zero  moment  and  displacement  at 
the  lag  hinge  where  a  hingeless  rotor  requires  zero  angular 
displacement  and  deflection  at  the  root  end.  In  both  cases, 
shear  and  moment  must  be  zero  at  the  blade  tip. 

Linear  beam  theory  was  utilized  as  a  means  of  verifying 
output  of  the  Myklestad  bending  program.  In  the  case  of  the 
flexbeam,  where  clamped  end  restoring  moments  were  needed, 

F7 

the  rotary  spring  constant  of  —  =  —  was  used  to  check  the 

Myklestad  output.  This  rotary  spring  constant  for  the  blade, 
namely  the  flexbeam  root-end  moment  per  unit  tip  angular 
deflection,  was  a  required  input  for  the  NPS  rotor 
simulation. 

Frequencies  of  oscillation  were  checked  against  Den 
Hartog's  text  "Mechanical  Vibrations"  [Ref.  Bl]  .  For 
example,  the  first  lag  frequency  of  the  hingeless  blade 
determined  by  the  Myklestad  program  was  checked  with  the 

cantilevered  beam  frequency  from  Den  Hartog,  a  =  3.52 
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The  plots  in  Figure  B3  shows  typical  output  of  the 

Myklestad  program.  The  third  subplot,  the  moment  as  a 
function  of  beam  longitudinal  position,  shows  a  comparison 
to  linear  theory.  The  example  beam  used  in  this  comparison 
has  a  length  of  10  inches,  a  weight  distribution  of  0.1 
lb/in,  and  a  constant  El  of  10000  lb*in2.  It  can  be  seen 

that  the  Myklestad  program  correctly  applies  the  four 
boundary  conditions  of  zero  slope  and  deflection  at  the 

root,  and  zero  shear  and  moment  at  the  tip.  It  should  be 
noted  that  a  small  correction  is  made  to  the  shear  at  the 
tip  of  the  beam.  This  is  done  due  to  the  final  Myklestad 

element  having  a  net  shear  at  the  tip  one  beam  element  away 
from  the  final  mass  element.  This  tip  correction  to  shear  is 
common  practice  in  applications  of  Myklestad' s  method. 
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Figure  B3  -  Myklestad  analysis  output 


The  following  table  compares  a  non-rotating 
cantilevered  beam  to  Den  Hartog's  closed  form  equations, 
verifying  the  frequencies  obtained  from  the  Myklestad 
program: 


Table  B1 


Linear  theory 

Myklestad 

Non-rotating  cantilevered  beam 

69.6484  Hz 

69.2986  Hz 

The  Myklestad  computer  code  follows  in  the  appendices. 
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APPENDIX  C  -  MYKLESTAD  METHOD  PROGRAM  IN  MATLAB 


function  myklestad 
clear 

format  compact 

%  Required  inputs: 

Ohm=0 ; 
tolb=0 . 1; 

(rad/sec) 

load  h66blade .mat  -ascii 
el=2 .5984 ; 

%  snubber  radial  location,  in 
radius=38 . 976; 

%  blade  radius,  in 

rb=fliplr (h66blade ( : , 1) ' *radius) ;  %  in 

wb=fliplr (h66blade ( : , 2 ) ' ) ;  %  lb  inA2 

EIb=fliplr(h66blade( : ,3) ' ) ;  %  lb/ in 

rf=rb (1:17) ; 
wf=wb{l : 17) ; 

EIf=EIb (1 : 17 ) ; 

figure(4) , plot (rf, Elf /1000) , grid 
xlabel ( 'Radius  (in) ' , 'FontSize' , 16) 
ylabel ( 'El ,  inA2Klb' , 'FontSize' , 16) 

title('EI  vs  Radial  Location  for  Flexbeam' FontSize' ,  16) 
muf= (wf/32 . 2) ; 
nfe=length(rf ) ; 

doutbd=[diff (rf) ,0] ; 
dinbd= [0, dif f (rf ) ] ; 

mEIf=sum( (dinbd+doutbd) /2 . *EIf ) / (rf (nfe) -rf (1)  )  ; 

disp(['  The  flexbeam  length  is:  ' , num2str (rf (nf e) -rf (1) ) , '  in']) 

disp ( [ '  The  mean  El  for  the  flexbeam  is:  ' , num2str (mEIf ) , '  lb*inA2']) 

%  obtain  beam  natural  frequency  and  root-end  transfer  matrix 
%  Finds  natural  frequency  and  dynamics  matrix  using 
%  Myklestad  method  for  beam  bending 

k=0  ; 

if  nfe~=length (Elf ) 

error ('The  El  vector  is  not  the  same  size  as  the  r  vector.') 
elseif  nfe~=length(muf ) 

error ('The  mass/length  vector  is  not  the  same  size  as  the  r 
vector. ' ) 
end 

omega_l=3 . 52*sqrt (mean(EIf ) / (mean(muf) * (rf (length(rf ) ) -rf (1) ) A4) ) ; 

%  rad/ sec 

omega= [0.9  1.1]*  omega_l ; 


%  Rotor  rad/ sec 

%  tolerance  in  acceptable  assumed  frequency 
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%******** *********WHILE  LOOP*** ****************** 
while  (abs (omega (1) -omega (2) ) >tolb) 
k=k+l; 

%  Th_tip  =  1,  all  others  zero 
tip= [1;  0;  0;  0]  ; 

state=mykloop (tip, omega (1) , rf ,muf , Elf , Ohm) ; 

%  compute  b  components  of  dynamics  matrix 
bTh=state (1,1) ;  bX-state (2 , 1 ) ;  bM=state ( 3 , 1 ) ;  bS=state ( 4 , 1 ) ; 
%b= [ bTh ,  bX,  bM] 

%  debug 

%plot state (2 ,r, state, omega (1) ) , pause 
%  debug 

%  X_tip  =  1,  all  others  zero 
tip— [0;  1/  0;  0]; 

state=mykloop (tip, omega (1) ,rf ,muf ,EIf ,Ohm) ; 

%  compute  a  components  of  dynamics  matrix 
aTh=state (1,1) ;  aX=state (2 , 1) ;  aM=state (3 , 1) ;  aS=state (4, 1) ; 
%a=[aTh,  aX,  aM] 

%  debug 

%plotstate (3 ,r, state, omega (1) ) , pause 
%  debug 

%  compute  determinant 

dynam= [bTh, aTh;bX, aX] ;  %  cantilevered  blade 


if  k==l 

omega=fliplr (omega) ; 
err=det (dynam) ; 
else 

err= [det (dynam) , err] ; 

%  use  linear  interpolation  for  next  guess  at  omega 

omega= [max (0, omega (1) -err (1) * ( (omega (2 ) -omega ( 1 ) ) / (err (2) - 
err (1) ) ) ) ,  omega] ; 

% [spline (err, omega, 0) ,  omega];  %  optional  spline 

method 

end 


end  %  while  loop 

%*****************TOILE  loop*  *  ******************* 


%  debugging 

disp ( [ '  ' ,num2str (k) , '  iterations  in  obtaining  root  transfer  matrix']) 

figure (3) , plot (omega (2 :k+l) , err , omega (2 : k+1) ,err, 'o') 

xlabel ( ' omega ' ) , ylabel ( ' determinant ' ) , grid, title ( ' Convergence  Plot ' ) 


Mj=-bM; 
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Kj=aM/rf (nfe) ; 
disp('  ' ) 

disp('  For  Simulation  Inputs,  the  Flexbeam  Produces') 

disp(['  the  following  values  for  Omega  =  ' ,num2str (Ohm*30/pi) , ' 

RPM' ] ) 

disp(['  Mj  =  ' ,num2str (Mj ) , '  in*lb/rad  deflection  at  the  joint']) 

disp(['  aS  =  ' ,num2str (aS) , '  So/Xt  for  root  transfer  matrix']) 

disp(['  Kj  =  ' ,num2str (Kj ) , '  lb/in  deflection  at  the  joint']) 

plots tat ( 2 , r f , state , omega ( 1 ) ) 
tip= [1;  0;  0?  0]; 

state=mykloop (tip, omega (1) , rf ,muf , Elf ,Ohm) ; 
bTh=state (1, 1) ;  bX=state (2 , 1) ;  bM=state (3 , 1) ; 
plots tat ( 3 , r f , state , omega ( 1 ) ) 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

load  h66cuff.mat  -ascii 

nel=length(h66blade( : ,1) ) ; 
rc=f liplr (h66cuf f ( : , 1) ' *radius) ; 
wc=fliplr (h66cuf f ( : , 2) ' ) ; 
ra=rb (17 : length (rb) ) ; 
wa=wb ( 17 : length (wb) ) ; 

rca=[rc,ra] ; 

%  in 

wca=  [wc, wa]  ; 

%  lb  inA2 

S=0 ; 1=0 ;M=0 ; 
for  n=2 : length (rca) 
dr=rca (n) -rca (n-1); 
wdr=mean(wca (n-1 :n) ) ; 
rcent=mean(rca(n-l:n) ) ; 

S=S+rcent*dr*wdr; 

I=I+rcentA2*dr*wdr; 

M=M+dr*wdr ; 
end 

Mb=-  (elA2*MA2-2*el*M*S+SA2)  /  ( -elA2*M+2*el*S-I)  ; 

Ms=  (SA2-I*M)  /  (-elA2*M+2*el*S-I)  ; 

R--  (el*S-I)  /  (-el*M+S)  ; 

disp(['  Mb  =  '  ,n\im2str  (Mb)  ,  '  slugs  for  blade  mass  at  position  R' ]  ) 
disp(['  Ms  =  ' ,num2str (Ms) , '  slugs  for  mass  at  snubber  position']) 
disp(['  R  =  ' , num2str (R) , '  inches  radial  location  for  mass  Mb']) 
disp ( '  ' ) 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

function  s tate=mykloop ( t ip , omega , r , mu , El , Ohm) 
%  next  step  subfunction  in  Myklestad  method 


%  in 

%  lb  inA2 
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%  delta-r  in 


nbe=length(r) ; 

1= [dif f (r) ] ; 1= [1  1 (length (1) )] ; 

V=1./EI; 

%  1/ (in*lb) 

U=(1.A2) . / (2*EI)  ;  % 

1/lb 

G= (1 . A3 )  . / ( 3  *EI ) ;  % 

in/lb 

T=fliplr  (cumsum(fliplr  (mu.*r)  )  )  *0hmA2;  %  lb 

%  tip  correction  to  shear: 

tip (4) =tip (4) + (mu (nbe-1 ) *1 (nbe-1) * (omegaA2+OhmA2) *tip(2) ) * .47; 
state= [zeros (4, nbe-1) , tip] ;  %  initialize  state 


for  n= (nbe-1) : -1 : 1 

%  theta 
%  deflection 
%  moment 
%  shear 


Th=state(l,n+1) ; 
X=state(2,n+1) ; 
M=state(3,n+1) ; 
S=state ( 4 , n+1 ) ; 


Thn  =  Th* (1+T(n) *U(n) ) -M*V(n) -S*U(n) ; 

Xn  =  X-Thn*l (n) +T (n) *Th*G (n) -M*U (n) -S*G (n) ; 

Mn  =  M+S*l  (n)  -T  (n)  *  (X-Xn)  ; 

Sn  =  S+mu(n)  *1  (n)  *  (omegaA2+OhmA2)  *Xn; 

%  [Thn;Xn;Mn;  Sn]  % 

debug 


state  ( :  ,  n)  =  [Thn;Xn;Mn;  Sn]  ; 

end 

%state 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function  [omegan, dynam, state] =natural (r ,mu, El, Ohm, articu, toler) 

%  Finds  natural  frequency  and  dynamics  matrix  using 

// 

%  Myklestad  method  for  beam  bending 

%  Robert  L.  King 
%  April  1998 

k=0  ; 

nbe=length(r) ; 
if  nbe-=length(EI) 

error (’The  El  vector  is  not  the  same  size  as  the  r  vector.') 
elseif  nbe~=length (mu) 
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error ('The  mass/length  vector  is  not  the  same  size  as  the  r 
vector. ' ) 
end 

if  articu==l  %  first  guess  at  1st  bending  freq  -  use  eqn  for 

uniform  beam: 

omega_l=15 . 4*sqrt  (mean (El)  /  (mean (mu)  *  (r  (length  (r) )  -r  (1) )  A4)  )  ; 
%  rad/sec 

else 

omega_l=3 . 52*sqrt  (mean (El)  /  (mean  (mu)  *  (r  (length  (r) )  -r  (1)  )  *4)  )  ; 
%  rad/ sec 

end 

omega= [0.9  1.1] *omega_l ; 

%*****************raiLE  L00P**** ******* ********** 
while  (abs (omega (1) -omega (2) ) >toler) 
k=k+l ; 

%  Th_tip  =  1,  all  others  zero 
tip= [1 ;  0;  0;  0]  ; 

state=mykll (tip, omega (1) , r,mu, El .Ohm) ; 

%  compute  b  components  of  dynamics  matrix 
bTh=state (1,1) ;  bX=state (2 , 1) ;  bM=state (3 , 1) ; 

%b= [bTh,  bX,  bM] 

%  debug 

%plotstate (2 ,r, state, omega (1) ) , pause 
%  debug 

%  X__tip  =  1,  all  others  zero 
tip= [ 0 ;  1 ;  0  ?  0] ; 

state=mykll  (tip,  omega  (1)  ,  r, mu, El, Ohm)  ; 

%  compute  a  components  of  dynamics  matrix 
aTh=state (1, 1) ;  aX=state ( 2 , 1 ) ;  aM=state ( 3 , 1 ) ; 

%a= [aTh,  aX,  aM] 

%  debug 

%plotstate (3 , r , state, omega (1) ) , pause 
%  debug 

%  compute  determinant 
if  articu==l 
dynam=  [bM,  aM;bX,  aX]  ; 
else 

dynam=  [bTh,  aTh;  bX,  aX]  ; 
end  %  if  statement 


%  articulated  blade 
%  cantilevered  blade 


if  k==l 

omega=fliplr  (omega)  ; 
err=det (dynam) ; 
else 

err= [det (dynam) , err] ; 

%  use  linear  interpolation  for  next  guess  at  omega 
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omega= [max (0, omega (1) -err (1) * ( (omega (2) -omega (1) ) / (err (2) 
err (1) )  )  )  ,  omega] ; 

% [spline (err, omega, 0) ,  omega];  %  optional  spline 

method 

end 


end  %  while  loop 

%******** *********raiLE  L00p*** ************ **★★★* 

disp(['  ' ,num2str (k) , '  iterations  in  natural .m.  '] ) 

%  debug 

omegan=omega  ( 1 )  ; 

%  debugging 

figure (1) , plot (omega (2 :k+l) , err , omega (2 : k+1 ) ,err, 'o' ) 

xlabel ( 'omega' ) ,ylabel { 'determinant' ) ,grid, title ( 'Convergence  Plot' ) 
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APPENDIX  D  -  MAPLE  WORKSHEET  FOR  COLEMAN  MODEL 


EQUATIONS  OF  MOTION  FOR  GROUND  RESONANCE  CONSIDERING  ONLY  INPLANE 
DEGREES  OF  FREEDOM 


DEFINE  COORDINATE  TRANSFORMATIONS 


Blade  undeformed  axis  --  Hub  : 

>  restart: 

>  with(linalg)  : 

Warning,  new  definition  for  norm 
Warning,  new  definition  for  trace 

>  psi :=Omega*t+Phi [k] ; 


psi  :=  Omega  t  +  Phi [k] 

>  Tl:=alpha->matrix(3,3, [1, 0, 0, 0, cos (alpha) , sin (alpha) ,0,- 
sin  (alpha)  , cos (alpha) ] ) ; 

T1  :=  alpha  ->  matrix (3,  3, 

[1,  0,  0,  0,  cos (alpha),  sin(alpha),  0,  -sin (alpha) ,  cos (alpha)]) 

>  T2:=alpha->matrix(3,3, [cos (alpha) , 0, - 
sin (alpha) , 0 , 1, 0 , sin (alpha) , 0 , cos (alpha)  ] )  ; 

T2  :=  alpha  ->  matrix (3,  3, 

[cos (alpha),  0,  -sin (alpha) ,  0,  1,  0,  sin (alpha) ,  0,  cos (alpha)]) 

>  T3 : =alpha->matrix (3,3, [cos (alpha) , sin (alpha) , 0 , - 
sin (alpha) ,cos (alpha) , 0, 0, 0,1] ) ; 

T3  :=  alpha  ->  matrix (3,  3, 

[cos (alpha),  sin(alpha),  0,  -sin(alpha),  cos (alpha),  0,  0,  0,  1]) 

>  diffl:=arg->map(diff ,arg, t) ; 

diffl  :=  arg  ->  map(diff,  arg,  t) 

>  Ml : =transpose (T3 (psi) ) ; 


[cos (%1) 

-sin(%l) 

0] 

[ 

] 

Ml  :=  [sin(%l) 

cos (%1) 

0] 

[ 

] 

[  o 

0 

1] 

%1  :=  Omega  t  +  Phi[k] 
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>  M2 :=transpose(T3 (zeta[k] (t) ) ) ; 


[cos (zeta[k] (t) ) 

-sin(zeta [k] ( t) ) 

0] 

[ 

] 

M2  :=  [sin(zeta[k] (t) ) 

cos (zeta [k] ( t) ) 

0] 

[ 

] 

[  0 

0 

1] 

Energy  of  rotor  blades 

Kinetic  energy  of  kth  rotor  blade  (TBk) 

>  rhoHI_I : =vector ( [u[l] (t) ,u[2] (t) , 0] ) ; 

rhoHI_I  :=  [u[l](t)/  u[2](t),  0] 

>  rhoBuH:=vector( [el,0,0] ) ; 

rhoBuH  :=  [el,  0,  0] 

>  rhoBuH_I : =multiply (Ml , rhoBuH) ; 

rhoBuH_I  :=  [cos(Omega  t  +  Phi[k])  el,  sin(Omega  t  +  Phi [k] )  el,  0] 

>  rhoPBd: =vector ( [R, 0 , 0 ] ) ; 

rhoPBd  :=  [R,  0,  0] 

>  rhoPBd_I:=multiply (Ml, M2, rhoPBd) ; 

rhoPBd_I  :=  [(cos(%l)  cos (zeta [k] ( t ) )  -  sin(%l)  sin ( zeta [k] ( t ) ) )  R, 
(sin(%l)  cos (zeta [k] (t) )  +  cos(%l)  sin(zeta[k] (t) ) )  R,  0] 

%1  :=  Omega  t  +  Phi[k] 

>  rho : =map ( simplify, matadd (rhoHI_I ,matadd (rhoBuH_I , rhoPBd_I) ) ) ; 
rho  :=  [u[l] (t)  +  cos(%l)  el  +  R  cos(%l)  cos (zeta [k] (t) ) 

-  R  sin(%l)  sin(zeta [k] (t) ) ,  u[2](t)  +  sin(%l)  el 
+  R  sin(%l)  cos (zeta [k] (t) )  +  R  cos(%l)  sin (zeta [k] (t) ) ,  0] 

%1  :=  Omega  t  +  Phi[k] 

>  V:=dif fl (rho) ; 

[/d  \ 

V  :=  [| —  u[l] (t) |  -  sin(%2)  Omega  el 
[\dt  / 
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-  R  sin(%2)  Omega  cos (zeta[k] (t) )  -  R  cos(%2)  sin(zeta[k] (t) )  %1 


-  R  cos (%2 )  Omega  sin(zeta[k] (t) )  -  R  sin(%2)  cos (zeta [k] (t) )  %1 

/d  \ 

,  | —  u[2] (t) |  +  cos{%2)  Omega  el 
\dt  / 

+  R  cos (%2 )  Omega  cos (zeta[k] (t) )  -  R  sin(%2)  sin(zeta [k] (t) )  %1 

-  R  sin(%2)  Omega  sin(zeta[k] (t) )  +  R  cos(%2)  cos (zeta [k] (t) )  %1 

] 

,  0] 

] 

d 

%1  :=  —  zeta[k] (t) 
dt 

%2  :=  Omega  t  +  Phi[k] 

>  TBk:  =l/2*mb[k]  *  (V[l]  A2+V[2]  A2)  ; 

///d  \ 

TBk  :=  1/2  mb[k]  | | | —  u[l] (t) |  -  sin(%2)  Omega  el 
\\\dt  / 

-  R  sin(%2)  Omega  cos (zeta[k] (t) )  -  R  cos(%2)  sin(zeta[k] (t) )  %1 

-  R  cos (%2 )  Omega  sin(zeta[k] (t) )  -  R  sin(%2)  cos (zeta [k] (t) )  %1 

\2  //d  \ 

|  +  | |--  u[2] (t) |  +  cos(%2)  Omega  el 

/  Wdt  / 

+  R  cos (%2)  Omega  cos (zeta [k] (t) )  -  R  sin(%2)  sin(zeta[k] (t) )  %1 

-  R  sin(%2)  Omega  sin(zeta[k] (t) )  +  R  cos(%2)  cos (zeta [k] (t) )  %1 

\2\ 

I  I 

/  / 

d 

%1  :=  —  zeta[k] (t) 
dt 

%2  :=  Omega  t  +  Phi[k] 


Potential  energy  for  kth  blade  (UBk) 
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>  UBkl : =l/2*Ke [k] *zeta [k] (t) A2 ;  #Linear  Elastic  Forces 


2 

UBkl  :=  1/2  Ke[k]  zeta[k](t) 

>  UBk2 :=l/4*Kd[k] *zeta [k] (t) A4;  #Duffing  Elastic  Forces 

4 

UBk2  :=  1/4  Kd[k]  zeta[k](t) 

>  UBk3  :  =l/4*Ks  [k]  *signum(zeta  [k]  (t) -z)  *  (zeta  [k]  (t)A2+zA2- 
2*zeta[k] ( t) *z) +l/4*Ks [k] *signum( zeta [k] (t) +z) * (-zeta[k] (t)A2-zA2- 
2*zeta [k] (t) *z) +l/2*Ks [k] *zeta[k] (t) A2+l/2*Ks [k] *zA2 ; 

UBk3  := 

1/4 

2  2 

Ks [k]  signum(zeta[k] (t)  -  z)  (zetafk] (t)  +  z  -  2  zeta[k] (t)  z) 

+  1/4 

2  2 

Ks [k]  signumfzeta [k] (t)  +  z)  (-zeta[k](t)  -  z  -  2  zeta[k] (t)  z) 

2  2 

+  1/2  Ks [k]  zeta [k] (t)  +  1/2  Ks[k]  z 

>  UBk : =UBkl+UBk2  +UBk3 ; 

>  UBk : =UBkl+UBk2 ; 

2  4 

UBk  :=  1/2  Ke [k]  zeta[k] (t)  +  1/4  Kd[k]  zeta[k] (t) 


Dissapative  function  for  kth  blade  (DBk) 


DBk: =l/2*Czeta [k] * (dif f (zeta [k] (t) , t) ) A2+Vzeta [k] * (diff (zeta [k] (t) ,t) ) A 
2*abs (dif f (zeta [k] ( t)  , t)  )  ; 

/d  \2 

DBk  :=  1/2  Czeta[k]  | —  zeta[k] (t) | 

\dt  / 

/d  \2  |  d  | 

+  Vzeta[k]  | —  zeta[k] (t) |  j  —  zeta[k] (t)  | 

\dt  /  j  dt  j 


Energy  of  hub 
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Kinetic  energy  of  hub  (TH)  /  Potential  energy  of  hub  (UH)  /  Dissapative 
function  of  hub  (DH) 


>  TF :=l/2*M[l]*(diff(u[l] (t) , t) ) A2+1/2*M [2 ] * (dif f (u [2 ] (t) ,t) ) A2; 

/d  \2  /d  \2 

TF  :=  1/2  M[l]  |  —  u [1]  (t)  |  +1/2M[2]  J  —  u [2]  (t)  | 

\dt  /  \dt  / 

>  UF : =1/2*K [1] *u [1] (t) A2+1/2*K[2] *u[2] (t)A2; 

2  2 
UF  :=  1/2  K[l]  u[l](t)  +  1/2  K[2]  u[2] (t) 


DF : =l/2*c [1] * (diff (u[l] (t) ,t) ) A2+l/2*c[2] * (diff (u[2] (t) ,t) ) A2+l/2*v[l] * 
(diff (u[l] (t) ,t) ) A2*abs(diff (u [1] (t) , t) ) +l/2*v[2] * (diff (u[2] (t) ,t) )A2*a 
bs (diff (u[2] (t) , t) ) ; 

/d  \2  /d  \2 

DF  :=  1/2  c[l]  |—  u[l]  (t)  |  +  1/2  c[2]  |—  u[2](t)| 

\dt  /  \dt  / 

/d  \2  |  d  | 

+  1/2  v [  1  ]  I—  u [  1  ]  ( t )  |  |  —  u[l](t)  j 

\dt  /  j  dt  j 

/d  \2  |  d  | 

+  1/2  v(2]  |—  u[2]  (t)  |  |  —  u[2]  (t)  j 

\dt  /  j  dt  j 


Generalized  forces  on  generalized  displacements 


>  F [1] :=0; 

F [1]  :=  0 

>  F[2] : =0 ; 

F [2]  :=  0 

>  F [3 ] : =u [ 1 ] ; 

F[3]  :=  u[l] 

>  F [4] : =u  [ 2 ]  ; 

F[4]  :=  u [2] 
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>  F [5] : =u [3 ] ; 


F  1 5  ]  u  [3  ] 


Derivation  of  equations  of  motion  using  Lagrange's  equation 


>  DOFF :  =  [u  [  1]  (t)  ,  U  [2]  (t)  ]  : 

>  DOFB:= [zeta [1] (t) , zeta [2] (t) ,zeta[3] (t) ] : 

>  DOF : = [ op ( DOFF ) , op ( DOFB ) ] : 

>  dDOF : =dif f 1 (DOF) : 

>  ddDOF : =dif f 1 (dDOF) : 

>  setA:={} :setB:={} :setC:={} : 

>  setD:={} :setE:={} :setF:={} : 

>  DOFq :  =  [ ] : dDOFq :  =  t ] : ddDOFq : = [ ] : 

>  for  i  from  1  to  vectdim(DOF)  do 

>  DOFq: = [op (DOFq) ,q[i] ] : 

>  dDOFq := [op (dDOFq) ,dq[i] ] : 

>  ddDOFq : = [ op ( ddDOFq) , ddq [ i ] ] : 

>  setA:=setA  union  {ddDOF [ i ] =ddDOFq [ i ]} : 

>  setB:=setB  union  {dDOF [ i] =dDOFq [i] } : 

>  setC:=setC  union  {DOF [i] =DOFq[ i] } : 

>  setD:=setD  union  {ddDOFq [i] =ddDOF [i] } : 

>  setE:=setE  union  {dDOFq [i] =dDOF [i] } : 

>  setF:=setF  union  {DOFqfi] =DOF [i] } : 

>  od: 

>  setl:=setA  union  setB  union  setC; 

2 
d 

setl  :=  { -  u[l] (t) 

2 

dt 

2 

d  d 

—  u[2](t)  =  ddq [  2  ]  ,  —  u[2] (t)  -  dq[2] ,  u[2] (t)  =  q[2], 

2  dt 

dt 

2 
d 

-  zeta [1] ( t) 

2 

dt 
2 

d  d 

-  zeta [2] (t)  =  ddq[4] ,  —  zeta[2](t)  =  dq[4] ,  zeta[2](t)  =  q[4] # 

2  dt 

dt 


d 

=  ddq [ 3 ] ,  —  zeta [1] ( t)  =  dq[3] ,  zeta[l](t)  =  q[3] , 
dt 


d 

=  ddq[l] ,  —  u[l] (t)  =  dq [ 1 ] ,  u[l](t)  =  q[l] , 
dt 
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2 

d  d 

-  zeta[3] (t)  =  ddq[5],  —  zeta[3] (t)  =  dq[5] ,  zeta[3] (t)  =  q[5] 

2  dt 

dt 

} 

>  set2:=setD  union  setE  union  setF; 

2 
d 

set2  :=  { ddq [ 1 ]  =  ---  u[l] (t) ,  dq[l] 

2 

dt 

2 

d  d 

ddq [ 2 ]  =  —  u[2] (t) ,  dq[2]  =  —  u[2] (t) ,  q[2]  =  u[2] (t) , 

2  dt 

dt 

2 
d 

ddq [ 3 ]  =  — 

2 

dt 

2 
d 

ddq [ 4 ]  =  - 

2 

dt 
2 

d  d 

ddq [ 5 ]  =  ---  zeta [3 ] ( t) ,  dq[5]  =  --  zeta[3]{t),  q [5]  =  zeta[3](t) 
2  dt 

dt 

} 

>  T:=TF: 

>  U : =UF : 

>  D1:=DF: 

>  for  i  from  1  to  vectdim(DOFB)  do 

>  T:=T+subs(k=i,TBk) : 

>  U:=U+subs (k=i,UBk) : 

>  Dl:=Dl+subs (k=i,DBk) : 

>  od: 

>  Temp:=subs (setl,T) : 

>  for  i  from  1  to  vectdim(DOF)  do 

>  tempi :=diff ( Temp , dDOFq [ i ] ) : 

>  temp2 : =subs ( set2 , tempi ) : 


d 

zeta[l](t),  dq[3 ]  =  --  zetaflHt),  q[3]  =  zeta[l](t) 

dt 


d 

zeta [2] (t) ,  dq [4]  =  —  zeta [2] (t) ,  q[4]  =  zeta[2] (t) 

dt 


=  —  u[l] (t) ,  q[l]  =  u[l] (t)# 
dt 
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>  temp3 :=diff (temp2,t) : 

>  LI : =subs ( setl , temp3 ) : 

>  L2 : =di f  f ( Temp , DOFq [ i ] ) : 

>  L3 :=diff (subs (setl, U) ( DOFq [i] ) : 

>  L4 : =dif f (subs (setl, Dl) ,dDOFq[i] ) : 

>  EOM [ i ] : =simplify (L1-L2+L3+L4-F [i] ) : 

>  od: 

>  sets : = (signum (1 , q [3 ] -z) =0 , signum(l , q [4] -z ) =0 , signum( 1 , q [5] - 

z) =0,signum(l,q[3] +z) =0, signum(l,q[4] +z) =0, signum (1 , q[5] +z) =0, abs (1, dq[ 
3 ] ) =0,abs (1, dq[4] ) =0 , abs (1, dq[5] ) =0, abs (1, dq[l] ) =0 , abs (1, dq[2] ) =0} : 

>  A:=matrix(vectdim(DOF) , vectdim(DOF) ) ; 

A  : =  array (1  ..  5,  1  ..  5,  []) 

>  for  i  from  1  to  vectdim(DOF)  do 

>  for  j  from  1  to  vectdim(DOF)  do 

>  A[i, j ] :=coef f (EOM[i] ,ddDOFq[j] ) : 

>  A[i, j ] :=subs (sets, A[i, j ] ) : 

>  od: 

>  od: 

>  Ax2dot:=multiply(A,ddDOFq) : 

>  f :=array(l. .vectdim(DOF) ) ; 

f  :=  array (1  . .  5,  [] ) 

>  for  i  from  1  to  vectdim(DOF)  do 

>  f [i] : =-simplify (EOM [i] -Ax2dot[i] ) : 

>  f [i] :=subs (sets, f [i] ) : 

>  od: 

>  xldot : = [ ] : xl : = [ ] : 

>  for  i  from  1  to  vectdim(DOF)  do  xldot := [op (xldot) ,x[i] ]  od: 

>  for  i  from  vectdim(DOF) +1  to  2*vectdim(DOF)  do  xl:=[op(xl) ,x[i] ]  od: 

>  setX:={} : 

>  for  i  from  1  to  vectdim(DOF)  do 

>  setX:=setX  union  {dDOFq[i] =xldot [i] } : 

>  setX:=setX  union  {DOFq[i] =xl [i] } : 

>  od: 

>  interface (labelling=false) ; 

>  A1 : =subs ( setX  ,op(A)); 

A1  :  = 

[M[l]  +  mb[l]  +  mb[2]  +  mb[3]  ,  0  , 

-mb[l]  R  cos (%1 )  sin(x[8] )  -  mb[l]  R  sin(%l)  cos(x[8])  , 

-mb [2]  R  cos (%2 )  sin(x[9])  -  mb[2]  R  sin(%2)  cos(x[9])  , 

-mb [3]  R  cos(%3)  sin(x[10])  -  mb[3]  R  sin(%3)  cos(x[10])] 

[0  ,  M[2]  +  mb[l]  +  mb[2]  +  mb[3]  , 

-mb[l]  R  sin(%l)  sin(x[8])  +  mb[l]  R  cos(%l)  cos(x[8])  , 
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-mb[2]  R  sin(%2)  sin(x[9] )  +  mb[2]  R  cos(%2)  cos(x[9])  , 

-mb [3]  R  sin(%3)  sin(x[10])  +  mb [3]  Rcos(%3)  cos(x[10])] 

[ 

[-mb[l]  R  cos(%l)  sin(x[8])  -  mb[l]  R  sin(%l)  cos(x[8])  , 

2 

-mb[l]  R  sin (%1)  sin(x[8])  +  mb[l]  R  cos(%l)  cos(x[8])  ,  mb[l]  R 

] 

,0,0] 

[ 

[-mb [.2]  R  cos(%2)  sin(x[9])  -  mb[2]  R  sin(%2)  cos(x[9])  , 

-mb [2]  R  sin(%2)  sin(x[9])  +  mb[2]  R  cos(%2)  cos(x[9])  ,  0  , 

2  ] 
mb[2]  R  ,  0] 

[ 

[ -mb [ 3 ]  R  cos (%3 )  sin(x[10])  -  mb[3]  R  sin(%3)  cos(x[10])  , 

-mb [3]  R  sin(%3)  sin(x[10])  +  mb[3]  Rcos(%3)  cos(x[10])  ,0,0 
2] 

,  mb [3 ]  R  ] 

%1  :=  Omega  t  +  Phi[l] 

%2  : =  Omega  t  +  Phi [ 2 ] 

%3  :=  Omega  t  +  Phi [3] 

>  fl:=subs (setX  , op(f)); 

[ 

fl  :=  t-K[l]  x[6]  -  c[l]  x[l]  -  v[l]  x[l]  |  x[l]  | 

2  2 

+  mb[l]  cos(%2)  Omega  el  +  mb[l]  R  cos(%2)  Omega  cos(x[8]) 

-  2  mb[l]  R  sin(%2)  Omega  sin(x[8])  x[3] 

2 

+  mb[l]  R  cos(%2)  cos(x[8])  x[3] 

2 

-  mb[l]  R  sin(%2)  Omega  sin(x[8]) 

+  2  mb[l]  R  cos(%2)  Omega  cos(x[8])  x[3] 

2  2 
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-  mb[l]  R  sin(%2)  sin(x[8])  x[3]  +  mb[2]  cos(%l)  Omega  el 

2 

+  mb [2]  R  cos(%l)  Omega  cos(x[9]) 

-  2  mb[2]  R  sin(%l)  Omega  sin(x[9])  x[4] 

2 

+  mb[2]  R  cos(%l)  cos(x[9])  x[4] 

2 

-  mb[2]  R  sin(%l)  Omega  sin(x[9]) 

+  2  mb[2]  R  cos(%l)  Omega  cos(x[9])  x[4] 

2  2 

-  mb[2]  R  sin(%l)  sin(x[9])  x[4]  +  mb[3]  cos(%3)  Omega  el 

2 

+  mb [3]  R  cos(%3)  Omega  cos(x[10]) 

-  2  mb[3]  R  sin(%3)  Omega  sin(x[10])  x[5] 

2 

+  nib [3]  R  cos(%3)  cos{x[10])  x[5] 

2 

-  mb[3]  R  sin(%3)  Omega  sin(x[10]) 

+  2  mb [3]  R  cos(%3)  Omega  cos(x[10])  x[5] 

2 

-  mb[3]  R  sin (%3 )  sin(x[10])  x[5]  ,  -c[2]  x[2]  -  K[2]  x[7] 

2 

-  v[2]  x[2]  |  x[2]  |  +  mb[2]  R  cos(%l)  Omega  sin(x[9]) 

+  2  mb [2]  R  sin(%l)  Omega  cos(x[9])  x[4] 

2  2 

+  mb [2]  R  cos ( % 1 )  sin(x[9])  x[4]  +  mb[3]  sin(%3)  Omega  el 

2 

+  mb[3]  R  sin(%3)  Omega  cos(x[10]) 

+  2  mb[3]  R  cos(%3)  Omega  sin{x[10] )  x[5] 

2 

+  mb[3]  R  sin(%3)  cos(x[10])  x[5] 

2 

+  mb [3]  R  cos(%3)  Omega  sin(x[10]) 

+  2  mb[3]  R  sin(%3)  Omega  cos(x[10])  x[5] 


132 


2  2 

+  mb[3]  R  cos(%3)  sin(x[10])  x[5]  +  mb[l]  sin(%2)  Omega  el 

2 

+  mb[l]  R  sin(%2)  Omega  cos(x[8]) 

+  2  mb[l]  R  cos(%2)  Omega  sin(x[8])  x[3] 

2 

+  mb [ 1]  R  sin (%2 )  cos(x[8])  x[3] 

2 

+  mb[l]  R  cos(%2)  Omega  sin(x[8]) 

+  2  mb [1]  R  sin(%2)  Omega  cos(x[8])  x[3] 

2  2 

+  mb[l]  R  cos(%2)  sin(x[8])  x[3]  +  mb [2]  sin(%l)  Omega  el 

2 

+  mb [2]  R  sin(%l)  Omega  cos(x[9]) 

2 

+  mb[2]  R  sin(%l)  cos(x[9])  x[4] 

+  2  mb[2]  R  cos  (%1)  Omega  sin(x[9])  x[4]  , 

2 

-2  Vzeta[l]  x[3]  |  x[3]  |  -  mb[l]  Omega  el  R  sin(x[8])  +  u[l] 

3  3 

-  Ke [1]  x [ 8]  -  Kd[l]  x[8]  -  Czeta[l]  x[3],  -Kd[2]  x[9] 

-  Czeta [2 ]  x[4]  -  Ke[2]  x[9]  -  2  Vzetaf2]  x[4]  |  x[4]  |  +  u[2] 

2 

-  mb [2]  Omega  el  R  sin(x[9])/  -2  Vzeta[3]  x[5]  |  x[5]  |  +  u[3] 

3 

-  Ke [3 ]  x [10]  -  Kd [3 ]  x[10]  -  Czeta[3]  x[5] 

2  ] 

-  mb[3]  Omega  el  R  sin(x[10])] 

%1  :=  Omega  t  +  Phi [2] 

%2  :=  Omega  t  +  Phi[l] 

%3  :=  Omega  t  +  Phi [3] 

>  realib{ fortran) : 

>  B:=augment (Al, fl) : 

>  f ortran(B, optimized); 
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fortran (B) 


APPENDIX  E  -  EQUATIONS  OF  MOTION  FOR  COLEMAN  MODEL 


mb[3]*R*sin(Omega*t+Phi[3])*sin(q[5])*dq[5]A2- 

mb[3]*R*sin(Omega*t+PW[3])*cos(q[5])*ddq[5]+mb[2]*R*sin(Omega*t+Pta[2])*Ome 

gaA2*sin(q[4])- 

2*mb[2]*R*cos(Omega*t+PM[2])*Omega*cos(q[4])*dq[4]+mb[2]*R*sin(Omega*t+Phi 
[2])*sin(q[4])*dq[4]A2-mb[2]*R*sin(Omega*t+Phi[2])*cos(q[4])*ddq[4]- 
mb[3]*cos(Omega*t+Phi[3])*OmegaA2*e  1  - 

mb[3]*R*cos(Omega*t+Phi[3])*OmegaA2*cos(q[5])+2*mb[3]*R*sin(Omega*t+Phi[3]) 

*Omega*sin(q[5])*dq[5]-mb[3]*R*cos(Omega*t+Phi[3])*cos(q[5])*dq[5]A2- 

mb[3]*R*cos(Omega*t+PW[3])*sin(q[5])*ddq[5]+mb[3]*R*sin(Omega*t+Phi[3])*Ome 

gaA2*sin(q[5])- 

mb[l]*R*cos(Omega*t+Phi[l])*OmegaA2*cos(q[3])+2*mb[l]*R*sin(Omega*t+Phi[l]) 

*Omega*sin(q[3])*dq[3]-mb[l]*R*cos(Omega*t+Phi[l])*cos(q[3])*dq[3]A2- 

mb[l]*R*cos(Omega*t+Phi[l])*sin(q[3])*ddq[3]+mb[l]*R*sin(Omega*t+Phi[l])*Ome 

gaA2*sin(q[3])- 

2*mb[l]*R*cos(Omega*t+Phi[l])*Omega*cos(q[3])*dq[3]+mb[l]*R*sin(Omega*t+Phi 

[l])*sin(q[3])*dq[3]A2-mb[l]*R*sin(Omega*t+Phi[l])*cos(q[3])*ddq[3]- 

mb[2]  *cos(Omega*t+Phi[2])*OmegaA2*e  1  - 

mb[2]*R*cos(Omega*t+Phi[2])*OmegaA2*cos(q[4])+2*mb[2]*R*sin(Omega*t+Phi[2]) 
*Omega*sin(q[4])*dq[4]-mb[2]*R*cos(Omega*t+Phi[2])*cos(q[4])*dq[4]A2- 
mb[2]*R*cos(Omega*t+Phi[2])*sin(q[4])*ddq[4]- 
mb[  1  ]  *cos(Omega*t+Phi  [  1  ]  )*OmegaA2*e  1  - 

2*mb[3]*R*cos(Omega*t+Phi[3])*Omega*cos(q[5])*dq[5]+K[l]*q[l]+v[l]*dq[l]*abs( 

dq[l])+M[l]*ddq[l]+l/2*v[l]*dq[l]A2*abs(l,dq[l])+c[l]*dq[l]+mb[l]*ddq[l]+mb[2]* 

ddq[l]+mb[3]*ddq[l]=0 

mb[l]*ddq[2]+mb[3]*ddq[2]-mb[3]*R*cos(Omega*t+Phi[3])*OmegaA2*sin(q[5])- 

2*mb[3]*R*sin(Omega*t+Phi[3])*Omega*cos(q[5])*dq[5]- 

mb[3]*R*cos(Omega*t+PM[3])*sin(q[5])*dq[5]A2+mb[3]*R*cos(Omega*t+Phi[3])*cos 

(q[5])*ddq[5]-mb[3]*R*sin(Omega*t+Phi[3])*OmegaA2*cos(q[5])- 

2*mb[3]*R*cos(Omega*t+Phi[3])*Omega*sin(q[5])*dq[5]- 

mb[3]*R*sin(Omega*t+Phi[3])*sin(q[5])*ddq[5]- 

2*mb[2]*R*sin(Omega*t+Phi[2])*Omega*cos(q[4])*dq[4]- 

mb[3]*R*sin(Omega*t+Phi[3])*cos(q[5])*dq[5]A2- 

mb[  1]  *sin(Omega*  t+Phi[  1  ])*OmegaA2*e  1  - 

mb[l]*R*sin(Omega*t+Phi[l])*OmegaA2*cos(q[3])- 

2*mb[l]*R*cos(Omega*t+Phi[l])*Omega*sin(q[3])*dq[3]- 

mb[l]*R*sin(Omega*t+Phi[l])*cos(q[3])*dq[3]A2- 

mb[l]*R*sin(Ornega*t+Phi[l])*sin(q[3])*ddq[3]- 

mb[l]*R*cos(Omega*t+Phi[l])*OmegaA2*sin(q[3])- 

2*mb[l]*R*sin(Ornega*t+Phi[l])*Omega*cos(q[3])*dq[3]- 

2*mb[2]*R*cos(Omega*t+Phi[2])*Omega*sin(q[4])*dq[4]- 
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mb[2]*R*sin(Omega*t+Phi[2])*cos(q[4])*dq[4]A2- 
mb  [2]  *R*sin(0mega*t+Phi[2]  )*  sin(q  [4] )  *ddq  [4]  - 
mb[2]*R*cos(Omega*t+Phi[2])*OmegaA2*sin(q[4])- 

mb[l]*R*cos(Omega*t+Phi[l])*sin(q[3])*dq[3]A2+mb[l]*R*cos(Omega*t+Phi[l]):,:cos 

(q[3])*ddq[3]-mb[2]*sin(Omega*t+Phi[2])*OmegaA2*el- 

mb[2]*R*sin(Omega*t+PM[2])*OmegaA2*cos(q[4])- 

mb[3]*sin(Omega*t+Phi[3])*OmegaA2*el- 

mb[2]*R*cos(Omega*t+Phi[2])*sin(q[4])*dq[4]A2+mb[2]*R*cos(Omega*t+Phi[2])*cos 

(q[4])*ddq[4]+M[2]*ddq[2]+mb[2]*ddq[2]+l/2*v[2]*dq[2]A2*abs(l,dq[2])+K[2]*q[2]+ 

c[2]*dq[2]+v[2]*dq[2]*abs(dq[2])=0 

2*Vzeta[l]*dq[3]*abs(dq[3])- 

u[l]+mb[l]*ddq[2]*R*cos(Omega*t+Phi[l])*cos(q[3])+mb[l]*RA2*ddq[3]- 

mb[l]*ddq[l]*R*cos(Omega*t+Phi[l])*sin(q[3])- 

mb[l]*ddq[l]*R*sin(Omega*t+Phi[l])*cos(q[3])- 

mb[l]*ddq[2]*R*sin(Omega*t+Phi[l])!|!sin(q[3])+mb[l]*OmegaA2*el*R*sin(q[3])+Vze 

ta[l]*dq[3]A2*abs(l,dq[3])+Ke[l]*q[3]+Kd[l]*q[3]A3+Czeta[l]*dq[3]=0 

mb  [2]  *  OmegaA2*e  1  *R*  sin(q[4]  )- 

mb[2]*ddq[l]*R*sin(Omega*t+Phi[2])*cos(q[4])+mb[2]*ddq[2]*R*cos(Omega*t+Phi[2 

])*cos(q[4])+mb[2]*RA2*ddq[4]- 

mb[2]*ddq[2]*R*sin(Omega*t+Phi[2])*sin(q[4])+Vzeta[2]*dq[4]A2*abs(l,dq[4])- 
u[2]+Ke[2]*q[4]+Kd[2]*q[4]A3+Czeta[2]*dq[4]+2*Vzeta[2]*dq[4]*abs(dq[4])- 
mb[2]  *ddq  [  1  ]  *R*cos(Omega*t+Phi[2]  )*  sin(q  [4]  )=0 

-u[3]-mb[3]*ddq[l]*R*cos(Omega*t+Phi[3])*sin(q[5])- 

mb[3]*ddq[l]*R*sin(Omega*t+Phi[3])*cos(q[5])+mb[3]*OmegaA2*el*R*sin(q[5])+mb 

[3]*RA2*ddq[5]- 

mb[3]*ddq[2]*R*sin(Omega*t+PW[3])*sin(q[5])+mb[3]*ddq[2]*R*cos(Omega*t+Phi[3 

])*cos(q[5])+2*Vzeta[3]*dq[5]*abs(dq[5])+Vzeta[3]*dq[5]A2*abs(l,dq[5])+Ke[3]*q[5]+ 

Kd[3]*q[5]A3+Czeta[3]*dq[5]=0 
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APPENDIX  F  -  SIMPLIFIED  ROTOR/ FUSELAGE  EQUATIONS  OF  MOTION 


Equations  with  viscous  and  hydraulic  damping  in  the  hub 
degrees  of  freedom,  x  and  y.  [Ref.  19] 


Mx  *x+Cx*x+Vx*x*\i^  +  Kx*x 
=  mb,  *R* (£'j  *sin('P1  +  £,)  +  (£!  +  £, )2  *cos('P,  +  £,) 

+  mb2  *R*(£2  *smC¥2+£2)  +  (Q  +  £2)2  *cos(%+C2) 
+  mb3  *R*(C 3  *sin('F3  +£'3)  +  (Q+£3)2  *cos('i'3  +^3) 

My  *y  +  Cy  *  j  +  Vy  *;y*|;y|  +  Xy  *  y 
=  mb,  *R* (£  * sin(%  +  £,)-(&  +  £, )2  * sinOF,  +  £ ) 

+  mb2  *R*(£2  *sin C¥2+£2)-(Cl+£2)2  *smQ¥2  +f2) 
+  mb3  *R*(C3  *sin(%  +f3)-(n  +  f3)2  *sin('¥3  +  C3) 


mbj  *R2  *£,  +C^l  *£,  +Ke ,  *£,  +mb ,  *Q2  * e * R * sin(^'1 ) 
=  mb,  *x*R* sin('P,  +  £,)  — mb,  *y*R* cos^  +£,) 


mb2  *R2  *^2  +C(2  *£2  +Ke2  *£2  +mb2  *Q2  *e*R*sin(£'2) 
=  mb2  *  3c  *  R  *  sinC'Pj  +  £2)  —  mb2  *  y*  R*  cos('P2  +£2) 


mb3  *R2  *g3  +C^3  *£3  +  Ke3  *£3  +mb3  *£22  *e*R*sin(£'3) 
=  mb3  * x * R * sin^j  +  £3)  —  mb3  *y*/?*cos('P3  +C3) 
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APPENDIX  G  -  5  BLADED  EQUATIONS  OF  MOTION  WITH  POLYNOMIAL 
SNUBBER  FUNCTIONS 

2 

mb[4]  R  sin(Omega  t  +  Phi [4])  Omega  sin{q[6]) 

-  2  mb [4]  R  cos (Omega  t  +  Phi [4])  Omega  cos(q[6])  dq[6] 

2 

+  mb [4]  R  sin(Omega  t  +  Phi [4])  sin(q[6])  dq[6] 

-  mb [4]  R  sin(Omega  t  +  Phi [4])  cos(q[6])  ddq[6] 

2 

-  mb [5]  cos (Omega  t  +  Phi [5])  Omega  el 

2 

-  mb[5]  R  cos (Omega  t  +  Phi [5])  Omega  cos(q[7]) 

+  2  mb [5]  R  sin (Omega  t  +  Phi [5])  Omega  sin(q[7])  dq[7] 

-  2  mb [5]  R  cos (Omega  t  +  Phi [5])  Omega  cos(q[7])  dq[7] 

2 

-  mb[2]  R  cos (Omega  t  +  Phi [2])  cos(q[4])  dq[4] 

-  mb[2]  R  cos (Omega  t  +  Phi [2])  sin(q[4])  ddq[4] 

2 

+  mb[2]  R  sin(Omega  t  +  Phi [2])  Omega  sin(q[4]) 

-  2  mb[2]  R  cos (Omega  t  +  Phi [2])  Omega  cos(q[4])  dq[4] 

2 

+  mb [2]  R  sin (Omega  t  +  Phi [2])  sin(q[4])  dq[4] 

-  mb[2]  R  sin(Omega  t  +  Phi[2])  cos(q[4])  ddq[4] 

2 

-  mb [3]  R  cos (Omega  t  +  Phi [3])  Omega  cos(q[5]) 

+  2  mb[3]  R  sin (Omega  t  +  Phi [3])  Omega  sin(q[5])  dq[5] 

2 

-  mb [3]  R  cos (Omega  t  +  Phi [3])  cos(q[5])  dq[5] 

-  mb[3]  R  cos (Omega  t  +  Phi [3])  sin(q[5] )  ddq[5] 

2 

+  mb [3]  R  sin (Omega  t  +  Phi [3])  Omega  sin(q[5]) 

-  2  mb[3]  R  cos (Omega  t  +  Phi [3])  Omega  cos(q[5])  dq[5] 
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2 

mb[l]  cos(Omega  t  +  Phi[l])  Omega  el 

2 

mb[l]  R  cos (Omega  t  +  Phi[l])  Omega  cos(q[3]) 

2  mb[l]  R  sin (Omega  t  +  Phi[l])  Omega  sin(q[3])  dq[3] 

2 

mb[l]  R  cos (Omega  t  +  Phi[l])  cos(q[3])  dq[3] 

mb [ 1]  R  cos (Omega  t  +  Phi[l])  sin(q(3])  ddq[3] 

2 

mb[l]  R  sin(Omega  t  +  Phi[l])  Omega  sin(q[3]) 

2  mb[l]  R  cos (Omega  t  +  Phi[l])  Omega  cos(q[3])  dq[3] 

2 

mb[l]  R  sin (Omega  t  +  Phi[l])  sin(q[3])  dq[3] 

mbtl]  R  sin (Omega  t  +  Phi[l])  cos(q[3])  ddq[3] 

2 

mb [2]  cos (Omega  t  +  Phi (2])  Omega  el 

2 

mb [2]  R  cos (Omega  t  +  Phi (2])  Omega  cos(q[4]) 

2 

mb [3]  cos (Omega  t  +  Phi [3])  Omega  el 
2  mb[2]  R  sin(Omega  t  +  Phi [2])  Omega  sin(q[4])  dq[4] 

2 

mb [3 ]  R  s in ( Omega  t  +  Phi [3])  sin(q[5])  dq[5] 

mb[5]  R  cos (Omega  t  +  Phi [5])  sin(q[7])  ddq[7] 

2 

mb[5]  R  sin(Omega  t  +  Phi [5])  Omega  sin(q[7]) 

2 

mb [5]  R  sin(Omega  t  +  Phi [5])  sin(q[7])  dq[7] 

mb [5]  R  sin (Omega  t  +  Phi [5])  cos(q[7])  ddq[7] 

mb [3]  R  sin (Omega  t  +  Phi [3])  cos(q[53)  ddq[5] 

2 

mb [4]  cos (Omega  t  +  Phi [4])  Omega  el 

2 

rab[4]  R  cos (Omega  t  +  Phi [4])  Omega  cos(q[6]) 
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+  2  mb[4]  R  sin(Omega  t  +  Phi [4])  Omega  sin(q[6])  dq[6] 

2 

-  mb[4]  R  cos (Omega  t  +  Phi[4])  cos(q[6])  dq[6]  -  u[l] 

+  mb[5]  ddqtl]  +  mb[4]  ddqtl] 

+  mb (2 ]  ddqtl]  +  M[l]  ddq[l]  +  mb[l]  ddqtl] 

+  v [ 1 ]  dq[l]  |  dq[l]  |  +  mb[3]  ddqtl]  +  K[l]  q[l]  +  c[l]  dq[l] 

-  mb[4]  R  cos (Omega  t  +  Phi[4])  sin(q[6])  ddq[6] 

2 

-  mb [5]  R  cos (Omega  t  +  Phi[5])  cos(q[7])  dq[7]  =  0 

v[2]  dq[2 ]  |  dq[2]  | 

-  2  mb [3]  R  cos (Omega  t  +  Phi [3])  Omega  sin(q[5])  dq[5] 

2 

-  mb [3 ]  R  sin(Omega  t  +  Phi[3])  cos(q[5])  dq[5] 

-  mb[3]  R  sin(Omega  t  +  Phi[3])  sin(q[5])  ddq[5] 

2 

-  mb[3]  R  cos (Omega  t  +  Phi [3])  Omega  sin(q[5]) 

-  2  mb[5]  R  cos (Omega  t  +  Phi [5])  Omega  sin(q[7])  dq[7] 

2 

-  mb[5]  R  sin(Omega  t  +  Phi[5])  cos(q[7])  dq[7] 

-  mb [5]  R  sin (Omega  t  +  Phi[5])  sin(q[7])  ddq[7] 

2 

-  mb [5]  R  cos (Omega  t  +  Phi [5])  Omega  sin(q[7]) 

-  2  mb[5]  R  sin (Omega  t  +  Phi [5])  Omega  cos(q[7])  dq[7] 

2 

-  mb [5]  R  cos (Omega  t  +  Phi[5])  sin(q[7])  dq[7] 

+  mb[5]  R  cos (Omega  t  +  Phi [5])  cos(q[7])  ddq[7] 

-  2  mb[3]  R  sin (Omega  t  +  Phi [3])  Omega  cos(q[5])  dq[5] 

2 

-  mb [3 ]  R  cos (Omega  t  +  Phi[3])  sin(q[5])  dq[5] 

+  mb[3]  R  cos (Omega  t  +  Phi [3])  cos(q[5])  ddq[5] 
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2 

-  mb[4]  sin(Omega  t  +  Phi [4])  Omega  el 

2 

-  mb[4]  R  sin(Omega  t  +  Phi [4])  Omega  cos(q[6]) 

-  2  mb [4]  R  cos (Omega  t  +  Phi [4])  Omega  sin(q[6] )  dq[6] 

2 

-  mb[4]  R  sin(Omega  t  +  Phi[4])  cos(q[6])  dq[6] 

-  mb[4]  R  sin(Omega  t  +  Phi[4])  sin(q[6])  ddq[6] 

2 

-  mb[4]  R  cos (Omega  t  +  Phi [4])  Omega  sin(q[6]) 

-  2  mb[4]  R  sin(Omega  t  +  Phi [4])  Omega  cos(q[6])  dq[6] 

2 

-  mb[4]  R  cos(Omega  t  +  Phi [4])  sin(q[6])  dq[6] 

+  mb[4]  R  cos (Omega  t  +  Phi [4])  cos(q[6])  ddq[6] 

2 

-  mb [5]  sin (Omega  t  +  Phi [5])  Omega  el 

2 

-  mb[l]  sin (Omega  t  +  Phi [ 1 ] )  Omega  el 

2 

-  mb[l]  R  sin (Omega  t  +  Phi[l])  Omega  cos(q[3]) 

-  2  mb[l]  R  cos (Omega  t  +  Phi[l])  Omega  sin(q[3])  dq[3] 

2 

-  mb[l]  R  sin(Omega  t  +  Phi [ 1 ] )  cos(q[3])  dq[3] 

-  mb[l]  R  sin(Omega  t  +  Phi[l])  sin(q[3])  ddq[3] 

2 

-  mb[l]  R  cos (Omega  t  +  Phi[l])  Omega  sin(q[3]) 

-  2  mb[l]  R  sin(Omega  t  +  Phi[l])  Omega  cos(q[3])  dq[3] 

2 

-  mb[l]  R  cos(Omega  t  +  Phi[l])  sin(q[3])  dq[3] 

+  mb[l]  R  cos (Omega  t  +  Phi[l])  cos(q[3])  ddq[3] 

2 

-  mb [2]  sin (Omega  t  +  Phi [2])  Omega  el 

2 

-  mb[2]  R  sin(Omega  t  +  Phi [2])  Omega  cos(q[4]) 
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-  2  mb[2]  R  cos (Omega  t  +  Phi [2])  Omega  sin(q[4])  dq[4] 

2 

-  mb[2]  R  sin(Omega  t  +  Phi[2])  cos(q[4])  dq[4] 

-  mb[2]  R  sin(Omega  t  +  Phi[2])  sin(q[4])  ddq[4] 

2 

-  mb[2]  R  cos (Omega  t  +  Phi [2])  Omega  sin(q[4]) 

-  2  mb [2]  R  sin (Omega  t  +  Phi [2])  Omega  cos(q[4])  dq[4] 

2 

-  mb[2]  R  cos (Omega  t  +  Phi [2])  sin(q[4])  dq[4] 

+  mb [2]  R  cos (Omega  t  +  Phi [2])  cos(q[4])  ddq[4] 

2 

-  mb [3]  sin (Omega  t  +  Phi [3])  Omega  el 

2 

-  mb  [5]  R  sin  (Omega  t  +  Phi  [5])  Omega  cos(q[7])  +  mb[2]  ddq[2] 

+  M[2]  ddq [ 2 ]  -  u[2]  +  K[2]  q[2] 

2 

+  c[2]  dq[2]  -  mb [3]  R  sin (Omega  t  +  Phi [3])  Omega  cos(q[5]) 

+  mb[l]  ddq [ 2 ]  +  mb[5]  ddq[2]  +  mb[3]  ddq[2]  +  mb[4]  ddq[2]  =  0 

mb [1]  ddq 1 2 ]  R  sin(Omega  t  +  Phi[l])  sin(q[3]) 

+  mb [1]  ddq [ 2 ]  R  cos (Omega  t  +  Phi [ 1 ] )  cos(q[3]) 

-  mb [1]  ddq [ 1 ]  R  cos (Omega  t  +  Phi[l])  sin(q[3]) 

-  mb[l]  ddq[l]  R  sin(Omega  t  +  Phi[l])  cos(q[3]) 

2  2 

+  mb[l]  Omega  el  R  sin(q[3])  +  mb[l]  R  ddq[3]  +  q [3 ]  Kpoly[5] 

+  dq [ 3 ]  Cpoly [5] 

3 

+  q [ 3 ]  Kpoly [2 ]  |  q[3]  |  +  q[3]  Kpoly[4]  |  q [3 ]  | 

3 

+  dq[3]  Cpoly [2 ]  |  q[3]  |  +  dq[3]  Cpoly[4]  |  q[3]  | 

2  4  5 

+  dq[3 ]  Cpoly [3 ]  q[3]  +  dq[3]  Cpoly [1]  q[3]  +  3  q[3]  Kpoly [1] 
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3 

+  2  q[3]  Kpoly[3]  -  u[3]  =  0 


2 

mb[2]  Omega  el  R  sin(q[4]) 

-  mb[2]  ddq[l]  R  cos (Omega  t  +  Phi (2])  sin(q(4] ) 

-  mb [2]  ddq[l]  R  sin (Omega  t  +  Phi [2])  cos(q[4]) 

-  mb[2]  ddq[2]  R  sin(Omega  t  +  Phi[2])  sin(q[4]) 

+  mb[2]  ddq[2]  R  cos (Omega  t  +  Phi [2])  cos(q[4]) 

3  5 

+  2  q[4]  Kpoly [3]  +  3  q[4]  Kpoly[l]  +  q(4]  Kpoly[5] 

2  4 

+  dq[4]  Cpoly [ 5 ]  +  mb[2]  R  ddq[4]  -  u[4]  +  dq[4]  Cpolyfl]  q[4] 

3 

+  q[4]  Kpoly [2 ]  |  q[4]  |  +  q[4]  Kpoly[4]  |  q[4]  | 

+  dq [4]  Cpoly [4]  |  q[4]  | 

3  2 

+  dq[4]  Cpoly [2]  |  q[4]  |  +  dq[4]  Cpoly[3]  q[4]  =  0 


2 

mb[3]  R  ddq [ 5 ] 

-  mb[3]  ddq [ 2 ]  R  sin(Omega  t  +  Phi [3])  sin(q[5]) 

+  mb[3]  ddq [ 2 ]  R  cos (Omega  t  +  Phi [3])  cos(q[5]) 

-  mb[3]  ddq[l]  R  cos (Omega  t  +  Phi [3])  sin(q[5] ) 

-  mb[3]  ddq[l]  R  sin(Omega  t  +  Phi [3])  cos(q[5]) 

2 

+  mb [3]  Omega  el  R  sin(q[5]) 

4  3 

+  dq[5]  Cpoly [ 1 ]  q[5]  -  u[5]  +  dq[5]  Cpoly[2]  |  q[5]  | 

+  dq[5]  Cpoly [4]  |  q[5]  |  +  q[5]  Kpoly[4]  |  q[5]  | 

3  2 

+  q [ 5]  Kpoly [2]  |  q[5]  |  +  dq[5]  Cpoly[3]  q[5] 

5  3 

+  3  q[5]  Kpoly [1]  +  2  q[5]  Kpoly[3]  +  q[5]  Kpoly[5] 
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+  dq[5]  Cpoly [5]  =  0 


3 

q[6]  Kpoly [ 2 ]  |  q[6]  | 

3 

+  q[6]  Kpoly [4]  |  q[6]  |  +  dq[6]  Cpoly[2]  |  q[6]  | 

2 

+  dq[6]  Cpoly [4]  |  q [6]  |  +  dq[6]  Cpoly[3]  q[6] 

4 

+  dq[6]  Cpoly [1]  q[6]  +  dq[6]  Cpoly[5]  +  q[6]  Kpoly[5]  -  u[6] 

5  3 

+  3  q[6]  Kpoly [1]  +  2  q[6]  Kpoly[3] 

+  mb [4]  ddq [ 2 ]  R  cos(Omega  t  +  Phi[4])  cos(q[6]) 

-  mb [4]  ddq [2]  R  sin (Omega  t  +  Phi [4])  sin(q[6] ) 

-  mb[4]  ddq [ 1 ]  R  sin(Omega  t  +  Phi[4])  cos(q[6]) 

-  mb[4]  ddq [ 1 ]  R  cos (Omega  t  +  Phi [4])  sin(q[6]) 

2  2 
+  mb [4]  Omega  el  R  sin(q[6])  +  mb[4]  R  ddq[6]  =  0 


2  2 
mb [5]  R  ddq [7]  +  mb[5]  Omega  el  R  sin(q[7]) 

-  mb[5]  ddq [ 1 ]  R  cos (Omega  t  +  Phi [5] )  sin(q[7] ) 

-  mb[5]  ddq[l]  R  sin(Omega  t  +  Phi [5])  cos(q[7]) 

-  mb[5]  ddq [ 2 ]  R  sin(Omega  t  +  Phi[5])  sin(q[7]) 

+  mb[5]  ddq [ 2 ]  R  cos (Omega  t  +  Phi[5])  cos(q[7]) 

2 

+  dq[7]  Cpoly [3]  q[7]  +  dq[7]  Cpoly[4]  |  q[7]  | 

3 

+  dq[7]  Cpoly [2]  |  q[7]  | 

3 

+  q[7]  Kpoly [4]  |  q [7]  |  +  q[7]  Kpoly[2]  |  q[7]  | 

4  5 

+  dq[7]  Cpoly [ 1 ]  q[7]  -  u[7]  +  3  q[7]  Kpoly[l] 
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3 

+  2  q[7] 


Kpoly [3]  +  q [7]  Kpoly[5]  +  dq[7]  Cpoly[5]  =  0 
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APPENDIX  H  -  5  BLADED  POLYNOMIAL  S-FUNCTION  INPUT  FILE 

function  [11,12,13,14,15,16] =poly5in 

%  This  m-file  serves  as  input  file  for  running  the  simulink 
%  S-function  simpleBp.m  for  the  5  bladed  RAH-66  Wind  Tunnel  Model. 

% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  Helo  Physical  and  Aerodynamic  parameters 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  Equivalent  Length  of  rotor  blade  (length) . 

R=26.5445;  %  inches 

%  Rotor  speed  (radians  per/ sec) 

Omega=870* (2*pi/60 ) ;  %  Run  # 

%  Hinge  offset  (length) . 

el=8 . 5235 ;  %  inches  (9.804,  7.6898  ??) 

%  Lead/ lag  stop  position  (rad) . 
z=pi; 

%  Mass  of  rotor  blades  (mass) . 

mb(l)=0. 013074;  %  slugs 

mb(2)=0. 013074; 

mb(3)=0. 013074; 

mb(4)=0. 013074; 

mb(5)=0. 013074; 

%  Effective  mass  of  fuselage. 

h=0.610;  %  hub  offset  in  meters 

M(l) =3 . 504/hA2*0 . 06852 ;  %  convert  fuselage  inertia  to  equivalent 
M (2 ) —1 . 495/h^2*0 . 06852 ;  %  mass  at  hub  &  convert  to  slugs 

%  Blade  azimuth  phase  angles  (radians) 

Phi (1) =0 ; 

Phi ( 2 ) =2  *pi / 5 ; 

Phi (3 ) =4*pi/5 ; 

Phi (4) =6*pi/5 ; 

Phi (5) =8*pi/5 ; 

%  Blade  Parameters 

%  Spring  stiffness  polynomial  for  lead-lag  (moment /radian) 

%  snubber  stiffness (lb/ in)  *  geogain ( in/rad)  *  (el- 

r_snub) (in) 
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gg=5 . 1121;  %  geometric  gain  of  snubber  movement  from  lag,  in/rad 
%  Kpoly= [0,0,0,0,79 . 71*gg+660 . 9141] * (el-2 .5984) ;  %  linearized  1995 

snubber 

Kbeam=[0  000  660.9141];  %  original  Myklestad  output 
%  Ksnub=  [38512436. 84*ggA5,  -  4491019. 58*ggA4, 183  33  6 . 70*ggA3  ,  - 

3147 . 25*ggA2 , . . . 

%  79.71*gg] *( el -2. 5984)  ; 

%  1995  snubber 

Ksnub= [336241084. 22*ggA5, -38666827. 45*ggA4, 1547897. 06*ggA3,- 
26483 . 84*ggA2 , . . . 

233. 32*gg] * (el-2. 5984) ; 

%  1992  snubber 
%  Kpoly=Kbeam+Ksnub; 

Kpoly=2 . 5* (Kbeam+Ksnub) ; 

%  Lead-lag  stop  spring  constants  (moment /radian) 

Ks (1) =0; 

Ks (2 ) =0 ; 

Ks (3 ) =0 ; 

Ks ( 4 ) =  0 ; 

Ks ( 5) =0 ; 

%  Linear  damping  in  lead-lag  (moment/ (rad/sec) ) 

%  snubber  damping(lb*s/in)  *  geogain (in/rad)  *  (el- 

r_snub) (in) 

%  Cpoly=[0,0,0,0,0.59*gg]* (el-2.5984) ;  % 

linearized  1995  snubber 

%  Cpoly= [116695. 68*ggA5, -14022. 78*ggA4, 609.0 *ggA3, - 

12 . 05*ggA2 , 0 . 59*gg] . . . 

%  *  (el-2 . 5984 ) ; 

%  1995  snubber 

Cpoly= [-842875 . 01*ggA5 , 87576 . 78*ggA4 , - 
2847 . 07*ggA3 , 25 . 01*ggA2 , 0 . 47*gg] . . . 

* (el-2 . 5984 ) ; 

%  1992  snubber 
Cpoly=2  *Cpoly ; 

%  Fuselage  Parameters 

%  Linear  springs  in  translation  (force/length) 

K (1) =1416 . 6/hA2* . 2248/39 . 37*12*1 . 14;  %  convert  rotary 

spring  to  lateral 

K (2 ) =1183 . 5/hA2* . 2248/39 . 37*12 ;  %  spring  at  hub  plane  & 

convert  to  lbf/in 


%  Linear  damping  in  translation  ( force/ ( length/sec) ) 

%  convert  rotary  damping  to  lateral 
%  damping  at  hub  plane  &  convert  to  lbf/in/s 
%  ?  c(l)=2*0.254/hA2*sqrt (18 .43/0 .254- (1.15*2*pi) A2)  *  . 06852; 
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%  ?  c(2)=2*3.75/hA2*sqrt(64. 05/3. 75- ( 0 . 65*2 *pi) ~2 ) * . 06852; 

%  c(l)=6. 58*0. 22481*39. 37/hA2; 

%  c(2)=10. 21*0. 22481*39. 37/hA2; 

c=M*2 . * [0 . 0467, 0 . 0542] . *sqrt (K. /M) ; 

%  Non-linear  damping  in  translation (force/ (length/sec) ^2) 

v(l)=0; 
v(2) =0; 

%  Initial  conditions 


%  Blade  lead- lag  displacement  (radians) 

xli=0 ; 
x2i=0 ; 
x3i=0 ; 
x4i=0; 
x5i=0; 


%  Blade  lead- lag  rates  (radians /sec) 

xrli=0 ; 
xr2i=0; 
xr3i=0 ; 
xr4i=0; 
xr5i=0; 

%  Fuselage  translational  displacements  (length) 

xXi=0.1; 
xYi=0 . 1 ; 

%  Fuselage  translational  rates  (length/sec) 

xrXi=0 ; 
xrYi=0; 

%  Create  input  arrays : 

11  =  [mb,  M] ; 

12  =  [R, Omega, el, z] ; 

13  =  Phi; 

14  =  [c,  v,  Cpoly] ; 

15  =  [Kpoly,  Ks,  K] ; 

16  =  [xrXi , xrYi , xrli , xr2 i , xr3 i , xr4i , xr5i , xXi , xYi , xli , x2 i , x3 i , x4i , x5i 
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APPENDIX  I  -  5  BLADED  POLYNOMIAL  S-FUNCTION 

function  [sys,  xO]  =  bladeSp ( t , x, u, flag, II, 12 , 13 , 14, 15, 16) 


% 

% 

% 

% 


% 

% 

% 


S- function  arguments: 


t 

x 

u 

flag 


time 

state  vector 
input  vector 

switch  used  by  numerical  integration  (simulation) 
routine  to  access  certain  parts  of  the  s- function 


S-function  input  parameters: 


[mb(l)  ,mb(2)  ,mb(3)  ,mb(4)  ,mb(5)  ,M(1)  ,M(2)  ] 

[R, Omega, el, z] 

[ Phi ( 1 ) , Phi ( 2 ) , Phi ( 3 ) , Phi (4) , Phi (5) ] ] 
tc(l)  ,c(2)  / v ( 1 )  , v(2 )  , 

Cpoly(l)  , Cpoly  (2)  , Cpoly (3)  /Cpoly(4)  , Cpoly (5)  ] 

[Kpoly (1)  , Kpoly (2 )  ,Kpoly(3)  ,Kpoly(4)  ,Kpoly(5)  , 

Ks ( 1 ) , Ks ( 2 ) , Ks ( 3 ) , Ks ( 4 )  , Ks ( 5 ) , 

K ( 1)  ,  K  (2 )  ] 

[xrXi , xrYi , xrli , xr2 i , xr3 i , xr4i , xr5i , 
xXi , xYi , xli , x2 i , x3 i , x4i , x5i ] 

%  S-function  to  represent  dynamics  of  5  bladed  coupled  rotor- 
%  fuselage  model  which  considers  only  inplane  degrees  of 
%  freedom,  i.e.,  x  and  y  translational  fuselage  degrees  of  freedom 
%  and  lead-lag  rotor  blade  degrees  of  freedom. 

% 

%  Explaination  of  variables: 

% - 

% 


% 

mb 

~> 

mass  of  blade 

% 

M 

-> 

effective  mass  of  fuselage 

% 

R 

-> 

distance  from  lead-lag  hinge  to  blade  center  of  mass 

% 

el 

-> 

blade  hinge  offset 

% 

Omega 

-> 

rotor  speed 

% 

z 

-> 

angle  at  which  blade  hits  stops 

% 

Phi 

-> 

blade  phase  angle  w.r.t.  azimuth  post ion 

% 

c 

-> 

fuselage  linear  damping 

% 

V 

-> 

fuselage  hydraulic  damping 

% 

Cpoly 

-> 

blade  damping  polynomial  coefficients 

% 

K 

~> 

effective  stiffness  of  fuselage  (landing 

gear  stiffness) 

% 

Kpoly 

-> 

blade  elastic  spring  constant  polynomial 

coefficients 

% 

Ks 

-> 

blade  stop  effective  spring  constant 

% 

% 

% 

% 


II 


12 


% 

% 

% 

% 

% 

% 

% 


13 


14 


15 


%  16 
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%  xr _ i  ->  initial  rate 

%  x _ i  ->  initial  displacement 

% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Define  input  parameters 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

mb=Il (1:5) ; 

M=I1 (6:7) ; 

R=I2 (1) ;0mega=I2 (2) ;el=I2 (3) ;z=I2 (4) ; 

Phi=I3 ; 

C=I4 (1 : 2 ) ; 

v=I4 (3:4) ; 

Cpoly=I4 (5:9) ; 

Kpoly=I5 (1:5) ; 

Ks=I5 (6:10) ; 

K=I5 (11:12)  ; 
xrXi=I6 (1) ;xrYi=I6 (2 ) ; 

xrli=I6 (3) ;xr2i=I6 (4) ;xr3i=I6(5) ;xr4i=I6(6) ;xr5i=I6(7) ; 
xXi=I6 (8) ;xYi=I6 (9) ; 

xli=I6 (10) ;x2i=I6 (11) ;x3i=I6(12) ;x4i=I6(13) ;x5i=I6 (14) ; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  S-function  flag  conditionals 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
if  flag  ==  0 

sys=[14, 0,14,7,0,0]  ; 

x0= [xrXi , xrYi , xrli , xr2i , xr3 i , xr4i , xr5i , xXi , xYi , xli , x2i , x3 i , x4i , x5i ] ; 

elseif  abs(flag)  ==  1 

t2  =  mb(l) *R; 

t3  =  Omega*t; 

t4  =  t3+Phi (1) ; 

t5  =  cos (t4) ; 

t6  =  sin(x(10) ) ; 

t7  =  t5*t6 ; 

t9  =  sin (t4) ; 

tlO  =  cos (x(10) ) ; 

til  =  t9*tl0; 

tl3  =  -t2*t7-t2*tll; 

tl4  =  mb (2 ) *R; 

tl5  =  t3+Phi (2 ) ; 

tl6  =  cos (tl5) ; 

tl7  =  sin (x(ll) ) ; 

tl8  =  tl6*tl7 ; 

t20  =  sin (tl5) ; 

t21  =  cos (x ( 11 ) ) ; 

t22  =  t20*t21 ; 

t24  =  -tl4*tl8-tl4*t22 ; 

t25  =  mb(3) *R; 

t26  =  t3+Phi (3 ) ; 
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t27  =  cos ( t26) ; 

t28  =  sin(x(12) ) ; 

t29  =  t27*t28 ; 

t31  =  sin(t26) ; 

t32  =  cos (x (12)  )  ; 

t33  =  t31*t32 ; 

t35  =  -t25*t29-t25*t33 ; 

t36  =  mb (4) *R; 

t37  =  t3+Phi (4) ; 

t38  =  cos ( 1 3 7 )  ; 

t39  =  sin (x (13 ) ) ; 

t40  =  t38*t39 ; 

t42  =  sin(t37) ; 

t43  =  cos (x(13) ) ; 

t44  =  t42  *t43 ; 

t46  =  -t36*t40~t36*t44 ; 

t47  =  mb (5) *R; 

t48  =  t3+Phi (5) ; 

t49  =  cos ( t48) ; 

t50  =  sin(x(14) ) ; 

t51  =  t49*t50 ; 

t53  =  sin(t48); 

t54  =  cos (x(14) ) ; 

t55  =  t53*t54; 

t57  =  -t47*t51~t47*t55 ; 

t62  =  0megaA2; 

t63  =  t27*t62 ; 

t66  =  t2  0* tl7 ; 

t67  =  x  (4)  A2  ; 

t70  =  t25*t31; 

t72  =  Omega*t28*x ( 5) ; 

t74  =  t42  *t39 ; 

t75  =  x ( 6 ) A2 ; 

t78  =  t42*t62 ; 

t81  =  t3  8*t62 ; 

t84  =  t38*t43 ; 

t87  =  t53  *t50 ; 

t88  =  x (7 ) A2 ; 

t91  =  t36*t38 ; 

t93  =  Omega*t43*x(6) ; 

t95  =  t36*t42 ; 

t97  =  Omega*t39*x(6) ; 

t99  =  t49  *t54 ; 

tl02  =  t53*t62 ; 

tl05  =  tl6*t62 ; 

tl08  =  t2  *t5 ; 

tllO  =  Omega*tlO*x (3 ) ; 

tll2  =  t9*t6 ; 

tll3  =  x (3 ) A2 ; 

tll6  =  t49*t62; 

tl20  =  t62  *el ; 

tl22  =  -v(l) *x(l) *abs (x(l) ) -K(l) *x(8) +t25*t63*t32-tl4*t66*t67 
2*t70*t72-t36*t74*t75-t36*t78*t39+t36*t81*t43+t36*t84*t75- 
t47*t87*t88+2*t91*t93-2*t95*t97+t47*t99*t88- 
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t47*tl02*t50+tl4*tl05*t21+2*tl08*tll0- 
t2*tll2*tll3+t47*tll6*t54+mb  (4) *t38*tl20; 


tl31 

= 

t47* t49 ; 

tl33 

= 

Omega*t54*x(7) ; 

tl35 

= 

t20*t62 ; 

tl38 

= 

t47*t53 ; 

t!40 

= 

Omega*t50*x(7) ; 

tl42 

tl6*t21 ; 

tl45 

= 

tl4*tl6 ; 

tl47 

= 

Omega*t21*x(4) ; 

tl49 

= 

tl4*t20 ; 

tl51 

= 

Omega*tl7*x(4) ; 

t!53 

= 

t9*t62; 

tl56 

= 

t5*t62 ; 

tl59 

= 

t31*t28; 

tl60 

= 

x  ( 5 )  A  2  ; 

tl63 

= 

t5*tl0 ; 

ties 

= 

t2*t9 ; 

tl68 

= 

Omega*t6*x(3) ; 

till 

= 

t27*t32 ; 

tl74 

= 

t31*t62; 

tl77 

= 

t25*t27 ; 

tl79 

_ 

Omega*t32*x(5) ; 

tl81 

= 

mb (5) *t49*tl20+mb(l) *t5*tl20+mb(2)*tl6*tl20+mb(3)*t27*tl20+2*tl31*tl33- 
tl4*tl35*tl7-2*tl38*tl40+tl4*tl42*t67+2*tl45*tl47+u ( 1) -2*tl49*tl51- 
t2*tl53  *t6+t2*tl56*tl0-t25*tl59*tl60+t2  *tl63  *tll3-2*tl66*tl68- 
c(l) *x(l) +t25*tl71*tl60-t25*tl74*t28+2*tl77*tl79; 
tl86  =  -t2*tll2+t2*tl63; 
tl89  =  tl4*tl42-tl4*t66; 
tl92  =  -t25*tl59+t25*tl71; 
tl95  =  -t36*t74+t36*t84; 
tl98  =  -t47*t87+t47*t99 ; 
t23  0  =  u(2)+tl4*tl8*t67- 

v(2) *x(2) *abs (x(2) ) +2*t91*t97+t36*t81*t39+2*t70*tl79+t36*t78*t43+t36*t4 
4*t75+2*tl77*t72- 

c (2) *x{2) +mb(2) *t20*tl20+t2*tl56*t6+t25*t63*t28+t25*t29*tl60+2*t95*t93+ 
t25*tl74*t32+t25*t33*tl60+t36*t40*t75+2*tl49*tl47 ; 
t265  = 

t47*t51*t88+t47*tll6*t50+2*tl45*tl51+t2*t7*tll3+2*tl38*tl33+t47*t55*t88 
+2*tl31*tl40+2*tl66*tll0+t2*tll*tll3+t47*tl02*t54+tl4*tl05*tl7+tl4*t22* 
t67+mb (1) *t9*tl20+tl4*tl35*t21+mb(4) *t42*tl20+mb(5) *t53*tl20+mb(3) *t31* 
tl20+2*tl08*tl68-K(2) *x(9) +t2*tl53*tl0 ; 
t267  =  R^2 ; 
t269  =  x(10)A2; 
t271  =  0; 

t274  =  abs (x(10)  )  ; 
t275  =  t274A2 ; 
t276  =  t275*t274; 
t287  =  t269^2; 
t299  =  el*R; 

t302  =  -t269*Kpoly (4 ) *t271/2-x ( 10 ) *Kpoly (2) *t276- 
x(10) *Kpoly(4) *t274-x(3) *Cpoly(2) *t276-x(3) *Cpoly (4) *t274- 
x(3) *Cpoly (3 ) *t269-x(3) *Cpoly (1) *t287+u (3) - 
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3 .EO/2 .E0*t269*Kpoly(2) *t275*t271-x (10) *Kpoly(5) -x(3) *Cpoly (5) - 
3*t287*x(10) *Kpoly (1) -2*t269*x(10) *Kpoly (3) -mb(l) *t62*t299*t6 ; 
t305  =  x(ll)A2; 
t306  =  t305A2 ; 
t318  =  abs (x (11) ) ; 
t319  =  t318A2 ; 
t320  =  0; 
t324  =  t319*t318 ; 

t336  =  u(4) -x(4) *Cpoly(l) *t306-3*t306*x(ll) *Kpoly (1) - 
2*t305*x(ll) *Kpoly (3) -x(ll) *Kpoly(5) -x(4) *Cpoly (5) -mb (2) *t62*t299*tl7- 
3 .EO/2 .E0*t305*Kpoly(2) *t319*t320-x (4) *Cpoly(2) *t324- 
x(ll) *Kpoly (2) *t324-x(ll) *Kpoly(4) *t318-t305*Kpoly (4) *t320/2- 
x(4) *Cpoly (3 ) *t305-x(4) *Cpoly (4) *t318; 
t339  =  x(12)A2; 
t340  =  t339A2; 
t345  =  abs(x(12)); 
t348  =  t345A2 ; 
t349  =  t348*t345 ; 
t356  =  0; 

t370  =  u(5)-x(5)*Cpoly(l)*t340-x(5)*Cpoly(3)*t339- 
x(5) *Cpoly (4) *t345-x(12) *Kpoly (2) *t349-x(12) *Kpoly (4) *t345- 
x(5) *Cpoly(2)*t349-t339*Kpoly(4)*t356/2-2*t339*x(12)*Kpoly(3)- 
3 .EO/2 . E0*t339*Kpoly (2 ) *t348*t356-mb(3) *t62*t299*t28-x(5) *Cpoly(5) - 
x(12) *Kpoly (5) -3*t340*x(12) *Kpoly(l)  ; 
t372  =  x(13)A2; 
t374  =  abs (x (13) ) ; 
t375  =  t374A2 ; 
t376  =  0; 
t381  =  t372A2; 
t390  =  t375*t374 ; 

t404  =  u (6) -3 .EO/2 .E0*t372*Kpoly (2) *t375*t376-x(13) *Kpoly(5) - 
x(6) *Cpoly (5 ) -3*t381*x ( 13 ) *Kpoly (1) -2*t372*x(13) *Kpoly (3 ) - 
mb(4) *t62*t299*t39-x(6) *Cpoly (2) *t390-x(6) *Cpoly(4) *t374- 
x(13) *Kpoly (4) *t374-t372*Kpoly (4) *t376/2-x ( 13 ) *Kpoly (2) *t390- 
x(6) *Cpoly (3) *t372-x(6) *Cpoly(l) *t381; 
t406  =  x(14)A2; 
t408  =  abs (x (14) ) ; 
t409  =  t408A2 ; 
t410  =  0; 
t421  =  t409*t408 ; 
t431  =  t406A2 ; 

t43  8  =  u(7) -3.E0/2 ,E0*t406*Kpoly(2) *t409*t410-x(7) *Cpoly(4) *t408- 
mb(5) *t62*t299*t50-t406*Kpoly (4) *t410/2-x(14) *Kpoly (2) *t421- 
x(14) *Kpoly (4) *t408-x(7) *Cpoly (2) *t421-x(7) *Cpoly (3 ) *t406- 
x (7 ) *Cpoly (5) -x(14) *Kpoly (5) -3*t431*x(14) *Kpoly(l)- 
2*t406*x ( 14) *Kpoly (3 ) -x (7) *Cpoly (1) *t431; 

B(l, 1)  =  mb ( 4 ) +mb ( 5 ) +mb ( 1 ) +mb ( 2 ) +mb ( 3 ) +M ( 1 ) ; 

B  (1 , 2 )  =  0; 

B (1 , 3 )  =  tl3 ; 

B (1 , 4)  =  t24 ; 

B(l, 5)  =  t35 ; 

B (1 , 6)  =  t46; 

B  (1 , 7 )  =  t57 ; 

B (1 , 8 )  =  tl22+tl81 ; 
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B  (2 , 1)  =  0; 

B (2 , 2}  =  mb ( 1 )  +mb ( 2 )  +mb ( 3 )  +itib ( 4 )  +mb ( 5 )  +M ( 2 )  ; 
B(2,3)  =  tl86 ; 

B (2 , 4)  =  tl89 ; 

B (2 , 5)  =  tl92  ; 

B (2 , 6)  =  tl95 ; 

B  (2 , 7 )  =  tl98 ; 

B (2 , 8)  =  t230+t265; 

B  (3 , 1)  =  tl3  ; 

B (3 , 2)  =  tl86; 

B  (3 , 3)  =  mb(l)*t267; 

B  (3 , 4)  =  0; 

B  (3 , 5)  =  0; 

B  (3 , 6)  =  0; 

B (3 , 7)  =  0; 

B  (3 , 8)  =  t302 ; 

B(4,l)  =  t24; 

B (4, 2)  =  tl89  ; 

B  (4 , 3 )  =  0; 

B ( 4 , 4)  =  mb(2) *t267; 

B (4, 5)  =  0; 

B  ( 4 , 6)  =  0; 

B  (4 , 7)  =  0; 

B  (4 , 8)  =  t336; 

B ( 5 , 1)  =  t35; 

B ( 5 , 2 )  =  tl92 ; 

B (5 , 3 )  =  0; 

B (5 , 4)  =  0; 

B(5, 5)  =  mb (3) *t267; 

B  (5 , 6)  =  0; 

B (5 , 7 )  =  0; 

B ( 5 , 8)  =  t370 ; 

B ( 6 , 1)  =  t46 ; 

B ( 6 , 2 )  =  tl95 ; 

B  (  6 , 3  )  =  0; 

B  ( 6 , 4 )  =  0; 

B(6,5)  =  0; 

B (6 , 6)  =  mb(4)*t267; 

B  ( 6 , 7)  =  0; 

B (6 , 8)  =  t404 ; 

B (7 , 1)  =  t57 ; 

B (7 , 2 )  =  tl98  ; 

B  (7 , 3 )  =  0; 

B (7, 4)  =  0; 

B (7 , 5)  =  0; 

B  (7 , 6)  =  0; 

B (7 , 7)  =  mb(5) *t267; 

B  (7 , 8)  =  1 4 3 8  ; 

%  Calculate  derivatives 

[m,  n]  =size  (B)  ; 

A1=B( : ,l:n-l) ; 
fl=B( : ,n) ; 
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sys=zeros (1,2  *m) ; 
sys (1:7) =Al\f 1 ; 
sys ( 8 : 14 ) =x ( 1 : 7 ) ; 

%  Output  states 

elseif  abs(flag)  ==  3 

sys (1:14) =x; 


else 


sys  =  [  ]  ; 


end 
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18.  Mr.  Ajay  Sehgal . 1 

1204  Wedglea  Dr. 

Bedford,  TX  76021 

19.  Dr.  JingYen . 1 

1812  Lakemont  Dr. 

Arlington,  TX  76013 

20.  Dr.  Richard  L.  Bennett . 1 

Bell  Helicopter  Textron,  Inc. 

P.O.  Box  482 
Fort  Worth,  TX  76101 
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21.  Mr.  Walter  Sonnebom . 1 

1936  Loma  Linda  Ct. 

Fort  Worth,  TX  761 12 

22.  Dr.  Andrew  Lemnios . 1 


Rensselaer  Polytechnic  Institute 
JEC  4010 
Troy,  NY  12180 

23.  Dr.  Daniel  Schrage . . . 1 

250  Drummen  Ct. 

Atlanta,  GA  30328 

24.  Dr.  C.  Eugene  Hammond . . 1 

100  Cove  Colony  Rd. 

Maitland,  FL  32751 

25.  Mr.  Robert  J.  Huston . 1 

105  Hex  Dr. 

Yorktown,  VA  23692-3012 

26.  Dr.  Osa  Fitch . 1 

39378  Sunnyside  Rd. 

Clements,  MD  20624 

27.  Mr.  Henry  Kelley . 1 

3092  N.  Riverside  Dr. 

Lanexa,  VA  23089-9403 

28.  Dr.  William  Warmbrodt . 1 

MS  T12-B 

NASA  Ames 
Moffett  Field,  CA  94035 

29.  Dr.  Wayne  Johnson . 1 

MS  T12-B 

NASA  Ames 
Moffett  Field,  CA  94035 


30.  LCDR  Christopher  Robinson . 1 

1318  Wakefield  Drive 

Virginia  Beach,  VA  23455 

31.  Mr.  Andrew  W.  Kerr . 1 

21950  McKean  Rd. 

San  Jose,  CA  95120 

32.  Mr.  Robert  Tomaine . 1 


Attn:  SFAE-AV-RAH-TV  (Tomaine) 
Bldg.  5681 
Redstone  Arsenal 
Huntsville,  AL  35898-5000 
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33.  Dr.  Richard  Carlson 
20999  Sarahills  Dr. 

Saratoga,  CA  95070 

34.  Dr.  Robert  L  Sierakowski 
2725  Lymington  Rd. 

Columbus,  OH  43220 

35.  Dr.  G.  A.  Pierce 
Georgia  Institute  of  Technology 
Dept,  of  Aerospace  Engineering 
Atlanta,  GA  30332-0150 

36.  Mr.  William  H  Weller 
11  Pine  Rd 

W  Hartford,  CT  06119 

37.  Dr.  Dewey  Hodges 
Georgia  Institute  of  Technology 
Dept,  of  Aerospace  Engineering 
Atlanta,  GA  30332-0150 

38.  Mr.  Robert  Blackwell 
103  OldZoarRd 
Monroe,  CT  06468 

39.  Mr.  Euan  Hooper 
729  Rhodes  Dr. 

Springfield,  PA  19064-1609 

40.  Mr.  Louis  Silverthom 
1336  North  Pasadena 
Mesa,  A Z  85201 

41.  Dr.  Friedrich  Straub 
1760  E.  Halifax  St. 

Mesa,  AZ  85203 

42.  Mr.  George  Price 
Sikorsky  Aircraft 
PO  Box  9729 
MS  S-330A 

Stratford,  CT  06497-9129 

43.  Dr.  R.  Lowey 
Georgia  Institute  of  Technology 
Dept,  of  Aerospace  Engineering 
Atlanta,  GA  30332-0150 

44.  Dr.  Farhan  Gandhi 
Penn  State/Aerospace  Eng. 

233  Hammond  Bldg. 


University  Park,  PA  16802 

45.  Dr.  Indeijit  Chopra . 1 

Aerospace  Engineering 

University  of  Maryland 
College  Park,  MD  20742 

46.  Dr.  Benson  Tongue . 1 

Department  of  Mechanical  Engineering 

6129  Etcheverry  Hall 
University  of  California 
Berkeley,  CA  94720-1740 
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